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three  dimensional,  discrete  model  of  the  human  spine,  torso 
and  head  was  developed  for  the  purpose  of  evaluating  mechanical 
response, in  pilot  ejection  and  it  was  developed  in  sufficient 
generality  to  be  applicable  to  other  body  response  problems,  such 
as  occupant  response  in  aircraft  crash  and  arbitrary  loads  on  the 
head-spine  system.  There  are  no  restrictions  on  the  distribution 
or  direction  of  applied  loads,  so  a wide  variety  of  situations  can 
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be  treated 


The  anatomy  is  modelled  by  a collection  of  rigid  bodies,  which 
represent  skeletal  segments  such  as  the  vertebrae,  pelvis,  head, 
and  ribs,  interconnected  by  deformable  elements,  which  represent 
ligaments,  cartilageneous  joints,  viscera,  and  connective  tissues. 

’’hniques  for  representing  other  aspects  of  the  ejection  environ- 
ment, such  as  harnesses  and  the  seat  geometry,  are  also  included. 

The  model  is  valid  for  large  displacements  of  the  spine  and  treats 
material  nonlinearities. 

The  basic  model  is  modular  in  format,  so  that  various  compon- 
ents may  be  omitted  or  replaced  by  simplified  representations. 

Thus,  while  the  complete  model  is  rather  complex  and  involves  sub- 
Iptantial  computational  effort,  various  simplified  models  are 
available  that  are  quite  effective  in  duplicating  the  response  of 
the  complete  model  within  a range  of  conditions.  Three  methods  of 
solution  are  available  for  the  analysis:  direct  integration  in  time 
I by  either  an  explicit,  central  difference  method  or  by  an  implicit, 
^trapezoidal  method,  and  a frequency  analysis  mdthod. 

j -^Results  are  presented  for  a variety  of  conditions,  such  as 
indifferent  rates  of  onset,  ejection  at  angles,  effects  of  lumbar 
purvature,  and  eccentric  head  loadings.  It  is  shown  that  large 
initial  curvatures  and  perfectly  vertical  acceleration  loadings  re- 
sult in  substantial  flexural  response  of  the  spine,  which  cause 
large  bending  moments.  It  is  further  shown  that  the  combination  of 
the  spine  s low  flexural  stiffness,  initial  curvature,  and  mass 
eccentricity  are  such  that  stability  cannot  be  maintained  in  a 10  g 
ejection  without  restraints  or  spine-torso-musculature  interaction^^ 

The  complete  models  were  used  mainly  to  study  the  effects  of  \ 
the  rib  cage  and  viscera  on  spinal  response.  The  flexural  stiffness 
of  the  torso  is  increased  substantially  by  a visceral  model,  even 
though  it  has  no  inherent  flexural  stiffness.  In  addition,  the 
viscera  provide  significant  reductions  in  the  axial  loads. 


Modal  analyses  were  performed  on  several  of  the  models  under 
various  conditions.  Numerous  flexural  natural  frequencies  under  10 
ops  were  found,  but  the  lowest  axial  frequency  is  of  the  order  of 
20  cps.  We  hypothesized  that  the  peaks  in  the  5-7  Hz  range  in 
iriving  point  impedances  observed  experimentally  in  axial  shaker 
table  measurements  result  from  pareunetric  excitation  of  the  flexural 
nodes. 
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CHAPTER  I 


INTRODUCTION 


1.  Objectives 


ii  I 


The  spine  is  the  primary  structural  element  for  transmitting 
forces  to  the  upper  torso  and  head  in  high  acceleration  environ- 
ments such  as  pilot  ejection.  Thus  in  the  study  of  ejection 
response,  it  is  common  to  model  the  element  for  force  transmission 
by  a bar  or  beam  and  to  neglect  the  torso  and  rib  cage.  These  bar- 
beam  models  have  evolved  into  two  general  classes:  the  so-called 
continuum  models,  in  which  the  bar  is  considered  as  homogeneous, 
and  the  discrete  models,  in  which  the  individual  vertebrae  are 
represented  as  rigid  bodies  and  are  connected  in  series  by  deform- 
able elements,  which  represent  the  intervertebral  disc  and  other 
connective  tissues.  These  two  types  of  models  are  in  fact  very 
similar  in  character,  for  if  the  scale  of  discretization  employed 
in  the  homogeneous  models  is  comparable  to  the  number  of  vertebral 
levels,  the  difference  equations  of  the  homogeneous  models  will  be 
very  similar  to  that  of  the  discrete  models.  The  primary  distinc- 
tion between  the  two  types  of  models  lies  in  the  possibility  of 
directly  using  disc  and  ligament  properties  in  the  discrete  models, 
whereas  the  continuum  models  require  determination  of  extrapolated 
material  properties,  which  represent  the  composite  behavior  of 
the  discs  and  vertebrae.  Both  the  discrete  models  and  homogeneous 
bar-beam  models  that  have  been  developed  so  far  have  been  restrict- 
ed to  one  or  two  dimensional  behavior. 

The  principal  objective  of  this  investigation  is  the  develop- 
ment of  a three  dimensional,  discrete  model  of  the  spine  and  head. 
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In  addition,  the  model  was  developed  in  a manner  so  that  other 
aspects  of  the  torso,  such  as  the  rib  cage  and  viscera,  could  be 
modelled  and  their  effects  on  the  behavior  of  the  spine  investi- 
gated. This  interaction  of  the  spine  with  the  torso  is  parti- 
cularly important  in  responses  which  involve  substantial  flexure 
of  the  spine,  for  the  flexural  stiffness  of  the  spine  is  very  low, 
and  as  shown  in  results  to  be  presented  subsequently,  are  not 
sufficient  to  insure  the  stability  of  the  spine  in  acceleration 
environments  commonly  found  in  pilot  ejection.  Significant  flexure 
may  be  induced  either  by  initial  curvatures  of  the  spine,  or  by 
asymmetric  properties,  such  as  asymmetric  mass  distribution-  Thus 
the  ability  to  investigate  the  behavior  of  the  spine  in  situations 
involving  substantial  bending  is  of  practical  importance. 

Because  the  .treatment  of  elements,  such  as  for  example,  the 
rib  cage,  in  sufficient  detail  to  accurately  represent  its  behav- 
ior in  a wide  variety  of  situations  involves  substantial  comput- 
ational effort,  the  model  has  been  developed  so  that  portions  of 
it  may  be  replaced  by  simplified  representations.  These  simpli- 
fied representations  are  quite  effective  in  a more  limited  range 
of  situations.  Thus,  the  rib  cage  can  be  replaced  by  an  equiva- 
lent beam  model,  and  a detailed  representation  of  the  cervical 
vertebrae  can  be  replaced  by  a single  beam  element.  These  simpli- 
fications provide  significant  savings  in  computer  time,  and  are 
therefore  quite  valuable  when  parametric  studies  are  undertaken. 

A general  description  of  the  characteristics  of  the  model  is 
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given  in  the  third  section  of  this  chapter.  The  details  of  the 
mathematical  formulation,  material  properties,  and  anatomical  repre 


sentation  are  given  in  the  next  three  chapters . Some  of  the  more  | 

interesting  and  significant  results  obtained  during  the  course 
of  this  study  are  then  described.  Finally,  the  input  data  formats 
for  both  the  dynamic  simulation  and  the  graphics  package  are  given 
in  the  Appendices. 

2.  Literature  Review  I 

I 

To  put  this  work  in  proper  prospective,  we  will  first  review  i 

i 

previous  models  of  the  spine,  using  the  customary  classification  | 

of  continuous  and  discrete  models.  Latham  (1957)  is  usually  , 

cited  as  the  first  to  develop  a mathematical  model  for  describing  ; 

the  dynamic  response  of  the  spine  to  +G  acceleration.  Latham's  | 

z y 

one  degree  of  freedom  model  consisted  of  the  rigid  masses  repre-  | 

senting  the  body  and  the  ejection  seat,  interconnected  by  a spring.  j 

It  was  developed  to  study  the  dynamic  overshoot  of  the  body  when 
seat  cushions  of  varying  resiliency  were  placed  between  the  pilot 
and  the  ejection  seat.  Also  included  in  Latham's  work  is  the 

t- 

first  study  concerning  the  natural  frequency  response  of  the  human  ' 

j, 

body  in  the  seated  position.  ^ 

Payne  (1961)  also  developed  a discrete,  one  degree  of  free- 

i 

dom  model  of  the  spine.  A rigid  mass  was  used  to  represent  the  [ 

head  and  upper  torso,  and  the  spine  was  modelled  as  a spring  with 
a dashpot  in  parallel.  The  stiffness  of  the  spring  was  chosen  to 
match  the  lowest  axial  natural  frequency  of  the  human  body  as 
predicted  from  the  lowest  peak  in  the  axial  driving  point  impedance 
measur  *ments.  Although  this  single  spring  model  could  not  predict 

9 


the  force  distribution  in  the  spine#  it  was  and  still  is  consid- 
ered an  accurate  representation  of  the  dyneimic  response  associated 
with  the  acceleration  profile  of  the  ejection  seat.  Subsequently, 
an  eight  (’‘^gree  of  freedom  model  was  developed  by  Toth  (1966)  . 

It  consisted  of  rigid  masses  representing  vertebrae  Til  through 
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L5  and  the  pelvis,  interconnected  by  springs  and  dampers  which 
represented  the  intervertebral  discs.  This  was  the  first  use  of 
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multiple  mass,  damped  spring  models  and  the  first  discrete  model 
to  idealize  individual  discs. 

Orne  and  Liu  (1970)  proposed  the  first  model  that  included 
the  shear  and  bending  resistance  of  the  intervertebral  disc.  The 
model  employed  a small  strain,  large  displacement  formulation. 


Each  of  the  vertebrae,  T1  through  L5,  was  represented  as  a rigid 
body  in  two  dimensional  space  with  three  degrees  of  freedom  per 
vertebra.  Spinal  curvature  and  variations  of  disc  stiffness  with 
vertebral  level  were  treated.  A three  parameter  viscoelastic 
force-deflection  relation  was  us>id  to  represent  the  material  pro- 
perties of  the  intervertebral  discs.  Orne  and  Liu  were  also  the 
firbL  Lo  model  the  iiieitidl  pioperties  by  assigning  to  each  motion 
segment  the  total  inertia  of  the  associated  segment  of  the  torso. 
Although  this  appears  somewhat  unreasonable  in  that  the  motion  of 
the  viscera,  because  of  its  low  shear  stiffness,  is  obviously 
quite  different  from  the  motion  of  the  vertebrae,  it  was  quite 
successful  in  duplicating  the  characteristics  of  experimentally  ob- 
served response  and  has  been  used  by  many  other  investigators. 

The  success  of  this  procedure  may  be  explained  in  terms  of  added 
masses  resulting  from  the  stiffness  of  the  viscera:  thus,  it  is 
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similar  to  the  "added  mass"  technique  used  to  analyze  the  vibra- 
tion of  structures  within  a fluid.  Also  included  in  the  model 
was  the  eccentricity  of  the  mass  center  for  each  motion  segment, 
which  was  assumed  to  be  uniform  along  the  spinal  column  with 
each  segment  having  the  same  inertial  properties.  The  model  did 
not  include  the  interactions  of  the  spine  with  the  torso, 
ejection  seat,  or  harness  apparatus.  Failure  to  represent  these 
interactions  in  a large  displacement  formulation  results  in  un- 
realistic deformed  configurations  of  the  spinal  column  and  may 
invalidate  the  force  distributions  predicted  by  the  model. 

Prasad  and  King  (1974)  extended  the  Orne  and  Liu  model  by 
including  the  articular  facet  interaction.  The  motivation  for 
this  extension  was  to  model  a secondary  load  path  in  the  spinal 
column  which  is  effected  by  the  articular  facets  as  indicated  by 
the  experimental  work  of  Prasad,  et  al.  (1973) . The  interaction 
of  the  articular  facets  was  modelled  by  two  springs,  one  limiting 
relative  rotations  and  the  other  limiting  the  relative  sliding  of 
adjacent  vertebrae. 

Stiffness  values  for  the  articular  facets  appear  to  have  been 
chosen  rather  arbitrarily,  since  no  reference  is  made  as  to  how 
the  axial  stiffness  was  determined  and  no  value  for  the  rotational 
stiffness  is  cited.  Of  the  axial  facet  stiffness  values  listed, 
the  largest  values  are  assigned  in  the  lumbar  region  and  are  of 
the  same  order  of  magnitude  as  the  disc  axial  stiffness.  Such 
large  stiffness  values  may  be  realistic  in  certain  directions, 
where  the  facet  effectively  imposes  a kinematic  constraint.  How- 
ever, the  deformation-resisting  character  of  the  facet  joint  in 
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other  directions  should  be  modelled  with  a much  lower  stiffness 
value,  as  pointed  out  by  Schultz,  et  al.  (1973) . Also  included 
in  the  model  is  an  auxiliary  force  representation  of  the  ejection 
seat  and  harness  interaction,  although  the  details  of  this  aspect 
of  the  model  were  not  described. 

A parallel  history  can  be  traced  in  the  homogeneous  (or  con- 
tinuum) models.  The  first  continuiun  model  was  proposed  by  Hess 
(1956),  who  included  only  axial  response.  Subsequently,  Moffat, 
et  al.  (1971)  included  both  axial  and  bending  response  by  using  a 
beam  type  model.  However,  the  analysis  was  restricted  to  small 
displacements . 

Recently,  Liu,  et  al.  (1973)  developed  bar-beam  models,  in- 
cluding large  displacements  in  the  analyses.  The  stiffness  pro- 
perties of  this  model  were  based  on  that  of  the  isolated,  liga- 
mentous spine  and  the  responses  they  exhibited  demonstrated  very 
large  deflections. 

3.  General  Description  of  Model 

The  model  represents  the  human  body  *^y  a collection  of  rigid 
bodies  interconnected  by  deformable  elements.  The  rigid  bodies 
are  used  for  the  modelling  of  bones,  while  the  defoinoaable 
elements  are  used  to  model  ligaments,  muscles  and  connective 
tissues.  The  treatment  of  bones  as  rigid  bodies  is  preferable 
from  both  the  viewpoint  of  numerics  and  modelling,  for  the  stiff- 
ness of  bones  is  usually  orders  of  magnitudes  greater  than  that 
of  connective  tissue,  so  that  if  both  are  modelled  as  deformable. 
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the  resulting  numerical  problem  is  poorly  conditioned.  However, 
long  slender  bones,  such  as  ribs,  may  be  modelled  as  deformable. 
The  deformable  elements  may  also  be  used  to  model  entities  exter- 
nal to  the  body,  such  as  restraint  systems  and  harnesses. 

For  purposes  of  describing  the  model,  it  is  worthwhile  to 
distinguish  between  the  following; 

1)  The  computer-based  method  of  solution,  or  mathematical  model, 
which  is  a rather  general  system  for  the  treatment  of  the  dyna- 
mics of  collections  of  rigid  bodies  interconnected  by  deformable 
elements , and 

2)  The  sife]t^ific  models  of  the  spine,  torso  and  ejection  system, 
which  constitute  a data  base  for  the  computer  system. 

We  will  first  describe  in  general  terms  the  mathematical 
model  employed  in  the  computer  simulation.  This  is  followed  by 
a general  description  of  the  data  sets  which  have  been  developed 
for  modelling  the  spine,  head,  and  torso  in  ejection  problems. 

4.  Mathematical  Model 

The  computer  procedure  is  basically  a matrix  structural 
analysis  technique,  which  serves  as  a versatile  framework  for 
constructing^ the  equations  of  motion.  The  progreim  enables  these 
equations  of  motion  to  be  integrated  in  time  by  either  explicit 
or  implicit  techniques,  or  analyzed  by  modal  procedures,  which 
give  the  natural  frequencies  and  modes  of  the  model.  The  formul- 
ation is  completely  three  dimensional  and  treats  arbitrarily 
lurge  rotations  and  displacements  of  the  rigid  bodies.  However, 


the  deformation  of  some  of  the  elements  is  restricted  to  be  moder- 
ately small.  Material  properties  may  be  linear  or  nonlinear  and 
linear  viscous  forces  can  be  included. 


Nodes  and  Coordinate  Systems.  Two  types  of  nodes  are  used: 

a)  primary  nodes,  each  of  which  has  six  degrees  of  freedom 
consisting  of  three  translations  and  three  rotations;  the  cen- 
troid of  a rigid  body  must  be  a primary  node; 

b)  secondary  nodes,  each  of  which  is  connected  through  a rigid 
body  to  a primary  node  and  which  thus  has  no  independent  degrees 
of  freedom. 

An  arbitrary  number  of  secondary  nodes  may  be  associated 
with  any  rigid  body,  and  they  serve  principally  as  a means  of 
connecting  deformable  elements  to  a rigid  body  at  a point  other 
than  the  centroid. 

The  configuration  of  the  model  is  described  by  the  position 
and  orientation  of  the  primary  nodes.  The  original  position  of 
node  I is  denoted  by  (i=l  to  3 representing  the  x,  y,  and  z 

components) ; the  new  position  are  obtained  by  adding  the  dis- 
placements u^j,  so 


^il  = 


ll  ll 


(1.1) 


The  orientation  of  a primary  node  is  described  by  a triad  of 
orthogonal  unit  vectors  ^21'  ^31'  rotate  with  the  node. 

In  order  to  describe  the  system,  we  will  define  three  types 
of  coordinate  systems; 

1.  a fixed,  global  set  of  coordinates  (x,  y,  z) , or  x^; 
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2.  body  coordinates  (x,  y,  z)^;  a set  of  body  coordinates  is 

associated  with  each  primary  node,  so  that  x,  y,  and  z coincide 
with  °2I'  ^31'  respectively  for  each  node.  The  origin 

of  the  x^j^  system  must  be  the  centroid  of  the  mass  at  node  I; 

3.  element  coordinates  (x,  y,  z) ; a set  of  element  coordinates 
is  associated  with  each  element,  and  the  element  coordinates 
rotate  and  translate  with  the  element  in  a manner  to  be  specified 
later.  The  x,  y,  and  z axes  are  associated  with  unit  vectors  e^^, 
62 , and  respectively  for  each  element. 

Model  Elements.  The  model  consists  of  the  following  elements: 

1.  rigid  bodies 

2.  spring  elements 

3.  beam  elements 

4.  hydrodynamic  elements 

5.  elastic  surfaces 


Rigid  Bodies.  Each  rigid  body  may  be  arbitrarily  oriented  in 
three  dimensional  space  and  may  undergo  arbitrarily  large  rota- 
tions and  translations.  The  centroid  of  the  rigid  body  is 
designated  a primary  node,  (see  Fig.  1) , its  coordinates  in 
space  define  the  position  of  the  rigid  body.  Each  rigid  body  has 
both  translational  and  rotational  inertia.  The  orientation  of 
the  rigid  body  is  described  by  the  triad  of  orthogonal  unit 
vectors  b^^,  b2,  and  b^.  The.  e vectors  must  coincide  with  the 
principal  axes  of  the  moment  of  inertia.  The  moments  of  inertia 
about  bj^,  b2,  and  b^  are  Ij^,  I2,  and  I^,  respectively.  In 
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Rigid  body  representation  and  coordinate  systems^ 
global  coordinates  (x,y , z) ; ^body  coordinates  (x,y,z) 
and  element  coordinates  (x,y,z). 


addition  to  the  primary  nodes,  any  number  of  secondary  nodes  may 
be  associated  with  the  rigid  body. 

Spring  Elements.  Spring  elements  are  deformable  elements  with 
only  axial  stiffness,  which  may  interconnect  any  two  nodes  of  the 
system.  A typical  spring  element  is  shown  in  Fig.  2. 


Figure  2.  Spring  Element 

The  element  may  be  connected  to  either  primary  or  secondary  nodes. 
The  axial  force  in  the  element  will  be  denoted  by  T,  with  T 
positive  in  tension;  the  elongation  is  designated  by  6.  The 
axial  force-elongation  law  is 

T = kj^6  + (1.2) 

where  either  k^  and  k_  may  be  zero.  A tension  cutoff 


<>< 


may  be  added  to  that  T = 0 whenever  6 < 0;  this  is  useful  for 
ligaments  and  other  elements  that  become  slack  whenever  the  elonga- 
tion is  negative. 


Beam  Element.  A beam  element  may  interconnect  any  two  nodes, 
which  may  be  either  secondary  or  primary  nodes.  Beam  elements 
include  axial  stiffness,  torsional  stiffness  and  bending  stiff- 


ness. The  resulting  nodal  forces  are  shown  in  Fig.  3:  f^j  and 

ak  a a a 

f , arise  from  axial  stiffness,  and  m,  _,  m _,  m _ and  m _ arise 
xJ  yi  zi  yj  zJ 

from  the  bending  stiffness  about  the  two  principal  coordinates  of 


the  cross-section,  y and  z,  and  m^^  arises  from  torsional  stiff- 
ness. For  all  moments,  the  right  hand  rule  sign  convention  is 
used  as  shown  in  Fig.  3. 


The  orientation  of  the  y and  z axes  is  given  by  including  a 
third  node  for  each  beam  element,  called  an  orientation  node. 


that  lies  in  the  y-x  plane  of  the  original  orientation  of  the  beam. 

There  are  two  available  modes  for  computing  the  forces  and 
moments  in  the  beeim.  In  the  first  mode,  functional  forms  are 
assumed  for  the  overall  response  of  the  beam;  these  are  not  con- 
sistent with  any  homogeneous  material  properties  but  allov.  the 
introduction  of  certain  nonlinearities.  The  forms  are: 


axial  force 

*XJ  = = 

optional  T=0if6<0  (1.3) 


M*<i 


bending  in  x-z  plane 
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(1.4b) 


For  linear  homogeneous  materials,  the  bending  constants  are  given 
through  standard  engineering  analysis  by 


E I 


= * l®2yl> 


(1.5) 


E I 


= 


where  E is  Young's  modulus,  il  the  length  of  the  element,  and  I the 
section  moduli,  which  are  respectively 


z'^dydz 


I 

zz 


u 


y^dydz 


(1.6) 


where  the  integral  is  over  the  cross-sectional  area  of  the  ele- 
ment, A.  The  shear  factor  is  given  by 
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(1.7) 


where  G is  the  shear  modulus  and  A the  effective  area  in  shear. 

s 

The  constants  a^  and  a.2  ure  included  to  permit  an  approximation 


to  a cubic  moment-curvature  behavior. 


The  second  method  of  computing  the  bending  moments  and 
axial  force  in  the  beam  requires  the  cross-section  of  the  beam  to 
be  defined  as  a thin-walled  member.  The  cross-section  is  defined 
by  the  coordinates  z^,  i = 1 to  NI,  and  the  shape  is  assumed 
to  be  prismatic,  so  that  y^,  are  constant  with  respect  to  x. 

If  this  mode  is  used,  the  moments  are  computed  directly  from  the 
axial  stresses  o by  numerical  integration.  An  arbitrary  stress 
strain  law  of  the  form 


a - k®e  + (1.8) 

may  be  used,  with  the  option  of  tension  cutoffs.  No  shear  cor- 
rections are  made.  This  mode  is  useful  for  modelling  elements 
such  as  the  walls  of  the  torso. 

In  both  modes,  the  torsional  resistance  is  taken  to  be  a 
linear  function  of  the  torsional  deformation  and  independent  of  the 
other  stresses  in  the  element,  i.e. 


As 


XlJ 


(i.y) 


. t 

where  m^  is  the  torque,  k the  torsional  spring  constant  and 
the  torsional  deformation.  The  shear  forces  are  always  obtainable 
from  equilibrium  so  no  force-deflection  law  is  necessary. 


Hydrodynamic  Element.  This  element  is  illustrated  in  Fig.  4.  The 
element  is  a pentahedron,  with  the  two  opposing  triangular  faces 
considered  to  be  rigid.  The  three  nodes  of  each  triangular  face 
must  therefore  be  associated  with  the  same  primary  node.  There 
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are  no  restrictions  on  the  geometry  of  the  element,  other  than 
that  the  initial  volume  of  the  element  be  positive:  this  is  in- 
sured by  numbering  the  nodes  appropriately. 

The  force  deflection  characteristics  of  this  element  are  ob- 
tained from  a linear  pressure-dilatation  relationship.  The 
pressures  are  transmitted  through  the  rigid  triangular  plates  to 
the  associated  primary  node  in  an  energetically  consistent 
fashion.  In  addition  to  the  linear  pressure-dilatation  stiffness, 
a linear  viscosity  is  available. 

This  element  is  useful  for  modelling  components  of  the  body 
that  exert  resistance  primarily  to  compressive  deformations. 
Because  of  the  presence  of  the  rigid  plates,  the  resistance  tends 
to  be  directed  through  a line  of  action  connecting  the  centroids 
of  the  two  triangular  surfaces.  Thus  it  is  useful  for  modelling 
articular  facets,  which  have  very  strong  directional  properties, 
and  the  viscera  that  effect  resistance  primarily  through  a verti- 
cal axis. 


Elastic  Planes.  An  assemblage  of  planes  may  be  prescribed  in  the 
model  to  represent  surfaces  of  the  pilot's  seat.  Each  plane  is 
described  by  locating  three  points  on  the  p.’.ane,  as  shown  in 
Fig.  5.  The  planes  restrain  the  motion  of  the  nodes  so  that  when 
a node  penetrates  the  plane,  a force  proportional  to  the  extent 
and  rate  of  penetration  is  applied  to  the  node  in  a direction  nor- 
mal to  the  plane. 

All  planes  are  considered  to  be  rigidly  linked  together.  The 
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motion  of  this  assemblage  of  planes  is  prescribed  through  either 
acceleration,  velocity,  or  displacement  histories. 


5.  General  Description  of  Models 

For  purposes  of  illustrating  how  the  mathematical  model  is 
used  to  represent  the  pilot's  anatomy,  we  will  here  describe  two 
representative  models  that  have  been  used  in  these  studies.  The 
first  model  is  restricted  to  the  isolated  throacolumbar  spine, 
the  cervical  spine,  the  head,  the  seatback  and  restraint  system. 

The  second  model,  in  addition  to  the  preceeding,  includes  a re- 
presentation of  the  rib  cage  and  viscera. 

The  first  model  is  graphically  depicted  in  Figs.  6 and  7, 
which  show  a back  view  and  a side  view  of  the  model  in  the  seated 
position,  respectively.  In  all  models  described  in  this  report, 
the  standard  orientation  for  the  global  coordinate  system  is  as 
follows:  the  z-axis  is  positive  vertically  upward,  the  y-axis  is 
positive  towards  the  back  and  the  x-axis  is  oriented  sideways; 
thus  the  y-z  plane  corresponds  to  the  sagittal  plane,  the  x-z 
plane  corresponds  to  the  frontal  plane  and  the  x-y  plane  corres- 
ponds to  the  horizontal  plane. 

The  graphical  depiction  in  Figs.  6 and  7 show  only  the  rigid 
bodies  employed  in  the  model.  Each  vertebra  and  the  head  is  re- 
presented as  a rigid  body.  The  configuration  of  these  rigid  bodies 
are  prescribed  by  the  initial  position  of  the  primary  nodes  in 
X,  y,  z space:  each  primary  node  must  coincide  with  the  mass  center 
of  the  rigid  body.  The  positions  of  the  primary  nodes  are 
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Figure  7 
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indicated  in  Fig.  7 by  plus  signs.  As  can  be  seen  from  the  figure, 
the  primary  nodes  in  many  cases  do  not  lie  within  the  vertebrae 
because  the  Liu,  et  al  (1973)  segment  data  was  used  to  represent 
the  inertial  properties  of  the  body  in  this  model,  so  that  each 
vertebra  is  associated  with  a segment  of  the  torso. 

The  conceptualization  used  in  modelling  the  inertial  properties 
of  the  human  body  differs  markedly  from  that  used  in  modelling  the 
stiffness  properties.  The  stiffness  model  considers  each  vertebra 
as  a rigid  body,  with  the  spring  elements  and  becun  elements  inter- 
connecting these  rigid  bodies  in  a manner  so  as  to  approximate 
force  deformation  characteristics  of  the  human  body.  On  the  other 
hand,  from  an  inertial  viewpoint,  each  rigid  body  represents  a seg- 
ment of  the  complete  torso,  and  each  vertebrae  is  considered  to  be 
rigidly  embedded  in  an  associated  segment  of  the  torso.  This 
corresponds  to  the  inertial  approximation  developed  by  Orne  and  Liu 
(1971)  and  used  by  Prasad  and  King  (1974) . The  more  complex  models 
do  not  use  this  approximation,  but  as  a consequence,  are  based  on 
less  reliable  data. 

In  the  thoracolumbar  spine,  each  pair  of  vertebrae  is  con- 
nected by  seven  spring  elements  and  one  beam  element.  The  inter- 
vertebral disc  is  represented  by  a beam  element,  which  joins  the 
geometrical  centers  of  the  endplates  of  each  pair  of  adjacent 
vertebrae.  The  spring  elements  represent  the  following  ligaments 
and  connective  tissues:  the  pair  of  spring  elements  which  connect 
the  transverse  process  tips  represent  the  intertransverse  liga- 
ments; one  spring  element,  which  connects  the  spinous  process  tips, 
represents  the  intra-  and  supra-spinous  ligaments;  a pair  of 


elements  which  connect  posterior  points  on  the  vertebral  bodies, 
represent  the  ligaments  flava;  two  spring  elements  are  used  to 
represent  the  articular  facets.  The  latter  are  short,  stiff  ele- 
ments and  are  primarily  intended  to  represent  the  kinematic  con- 
straints resulting  from  the  facets.  All  of  these  elements  inter- 
connect secondary  nodes  on  the  rigid  body.  In  addition,  the  pri- 
mary nodes  are  connected  by  additional  beam  elements  which  repre- 
sent the  stiffness  of  the  torso  and  rib  cage.  Cubic  moment 
curvature  relations  are  used  in  these  elements  so  that  their 
bearing  on  small-displacement  response  is  neglible.  Details  as  to 
the  locations  of  the  nodal  points  and  the  material  properties  of 
the  deformable  elements  may  be  found  in  Chapter  IV. 

In  the  cervical  spine,  adjacent  vertebrae  are  connected  only 
by  elements  representing  the  disc,  the  interspinous  ligaments, 
and  the  articular  facets.  The  discs  are  represented  by  bectm  ele- 
ments, the  ligaments  by  spring  elements,  the  articular  facets  are 
represented  by  hydrodynamic  elements.  The  triangular  endplates 
of  the  hydrodynamic  element  may  be  seen  in  Figs.  6 and  7.  Because 
these  elements  have  resistance  primarily  through  a line  joining 
the  centroids  of  the  two  opposing  triangular  facets,  these  ele- 
ments are  more  effective  in  representing  the  directional  proper- 
ties of  articular  facets  than  spring  elements.  The  use  of  these 
elements  for  the  representation  of  the  facets  would  also  be 
desirable  in  the  lumbar  and  thoracic  spines,  but  the  procurement 
of  data  for  the  location  of  the  facet  planes  in  these  portions  of 
the  spine  has  not  been  completed. 

The  head  is  a single  rigid  body  joined  to  C2  by  a beam 
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element  (Cl  was  not  included  in  the  model) . In  all  the  simulations 
studied  here,  the  helmet  was  assumed  to  move  exactly  like  the  head, 
so  that  if  a helmet  was  included  in  the  study,  the  moment  of  in- 
ertia and  mass  of  the  helmet  was  simply  added  to  that  of  the  head. 

The  seatback  in  this  model  is  a plane  surface,  which  is 
vertically  aligned  and  the  bottom  of  the  seat  is  horizontal.  The 
seat  constrains  the  motion  of  the  rigid  bodies  only  when  they  come 
in  contact.  The  definition  of  the  seatback  is  quite  arbitrary,  as 
long  as  it  can  be  described  as  a series  of  planes,  so  that  alterna- 
tive seatback  designs  can  be  studied  by  the  model  by  simply 
altering  the  description  of  these  planes. 

The  restraint  system  in  this  model  consists  of  4 springs,  3 
connecting  the  vertebrae  Tl,  T2,  and  T3  with  the  seatback,  the 
other  connecting  a secondary  node  on  the  pelvis  with  another  point 
on  the  seatback.  The  upper  restraint  belt  is  represented  by  three 
springs  to  reduce  the  shear  deformation.  Again  the  positions  of 
these  nodes  are  indicated  in  Fig.  7.  The  orientations  and  method 
of  interconnection  for  these  elements  is  completely  arbitrary  so 
that  other  harness  systems  can  be  modelled.  However,  important 
aspects  such  as  friction  and  the  actual  details  of  the  geometry  of 
the  restraint  system  have  not  yet  been  included. 

One  of  the  more  complex  models  is  represented  in  Figs.  8 and 
9,  which  show  the  back  and  side  views,  respectively.  The  major 
aims  of  this  model  are  the  separation  of  the  inertial  aspects  of 
the  torso  from  that  of  the  spine,  the  inclusion  of  certain 
structural  aspects  of  the  rib  cage,  and  the  addition  of  an  inde- 
pendent load  path  through  the  viscera.  These  aims  were  implemented 


Figure  9.  Side  view  of  spine-torso  model  with  rib  cage. 
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as  follows. 


The  rib  cage  is  represented  by  a system  of  rigid  bodies  and 
deformable  elemevts  which  includes  separate  rigid  bodies  for  each 
of  the  ribs  and  the  sternum.  Since  each  rib  is  modelled  as  a 
rigid  body,  the  deformation  of  the  thorax  as  a whole  results  from 
the  rotation  of  the  ribs  and  the  deformation  of  the  costo-sternal 
cartilage.  Each  rib  is  connected  to  two  vertebrae  by  means  of 
three  deformable  elements,  which  represent  the  costo-vertebral 
joint.  These  deformable  elements  have  been  chosen  so  that  the 
directional  properties  of  the  joint  are  represented  and  an  axis 
of  great  rotational  flexibility  was  included.  The  ribs  are  con- 
nected to  the  sternum  through  the  costo-sternal  joint  by  a deform- 
able element,  which  represents  the  deformability  of  the  costo- 
cartilage.  This  model  is  thus  quite  adequate  for  representing  the 
additional  bending  stiffness  of  the  torso  that  is  provided  by  the 
rib  cage;  on  the  other  hand,  it  is  not  suitable  for  representing 
frontal  impact  where  significant  deformations  of  the  rib  itself 
may  take  place.  For  the  latter,  it  would  be  necessary  to  represent 
the  deformation  of  the  neck  of  the  rib  by  modelling  it  by  a beam 
element. 

The  viscera  are  represented  by  a stack  of  hydrodynamic 
elements,  which  are  illustrated  in  Fig.  1C.  The  hydrodynamic 
elements  have  stiffness  only  when  deformed  axially,  so  that  this 
column  does  not  have  any  resist.*  to  shear.  However  because 
rigid  endplates  are  included  between  each  vertical  layer  of  hydro- 
dynamic  elements,  the  system  does  resist  bending  and  maintains 
coherence  in  response  to  transverse  loads.  The  very  bottom  plate 
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of  the  viscera  is  connected  to  the  pelvis,  whereas  the  uppermost 
plate  is  connected  to  ribs  10  on  the  left  and  right  side;  this 
interconnection  represents  the  transfer  f the  load  to  the 
diaphgram.  No  hydrodynamic  elements  are  included  within  the  thorax. 

The  inertial  properties  of  this  model  were  obtained  by  sub- 
dividing the  mass  of  each  segment  of  the  torso  between  the  spine 
and  the  ribs  and  sternum  in  the  thoracic  regions,  and  between  the 
spine  and  the  viscera  in  the  lumbar  region.  The  distribution  in 
each  segment  was  chosen  so  that  the  total  mass  of  each  segment 
corresponds  to  the  data  of  Liu,  et  al,  and  so  tl  t the  moment  of 
inertia  of  the  masses  of  the  components  in  each  segment  have  a 
moment  of  inertia  equal  to  that  measured  by  Liu,  et  al.  Because 
the  mass  of  each  body  segment  is  partitioned  into  the  inertia 
associated  with  the  spine  and  the  inertia  of  the  thorax,  the 
rotation  of  a body  segment  may  differ  from  the  rotation  of  an 
embedded  vertebral  body.  The  model  of  the  thoracolumnar,  cervi- 
cal spine,  head,  seatback,  restraints  and  pelvis  are  identical 
to  that  of  the  previous  model.  Both  the  details  of  the  geometry 
and  material  properties  may  be  found  in  Chapter  IV. 
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MATHEMATICAL  FORMULATION  OF  MODEL 
1 . Nomenclature  and  Coordinate  Transformations 


A general  description  of  the  modeling  techniques  has 
been  given  in  Chapter  1.  In  this  chapter,  the  detailed  equations 
and  mathematical  procedures  for  the  formulation  and  solution  of  the 
governing  equations  will  be  presented. 

The  coordinate  systems  have  already  been  described  in  Chapter  1. 
In  addition  to  a global  coordinate  system  (x,y,z),  local  coordinates 
for  each  element,  (x,y,z)„,  and  for  each  node,  (x,y,z)_,  are  used. 

The  unit  vectors  for  these  coordinate  systems  are  (e, ,e_,e_)„,  for 
the  coordinate  system  of  element  E and  ^ for  the  coordinate 

system  of  rigid  body  I . 

The  unit  vectors  and  immediately  define  the  rotational 
transformation  of  any  vector  components  between  the  coordinate 
systems.  Thus,  if  we  consider  a vector  A with  global  components 
(A^,A^,A^) , body  coordinate  components  ,A  ) and  element  coordin- 


(A 

X 

,A  ,A 
y z 

) , we 

have 

b, 

lx 

"2x 

^3x 

= 

ly 

^2y 

Iz 

"2z 

'’3Z. 

= tX]{A} 


(2.1) 
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where  b.  , b.  , b.  are  the  global  components  of  the  body  vectors. 


Similarly 


1 i 
li 


®lx 

®2x 

®3x 

®iy 

®2y 

®3y 

®lz 

®2z 

3z 

[vHa} 


(2.2) 


where  ®iy'  ®iz  Slobal  components  of  the  element  vectors. 

Also 


{A}  = 

[>]^{A} 

(2.3) 

A. 

{A}  = 

[yl'{A} 

(2.4) 

The  translational  motion  of  the  system  is  described  by  the  displace- 
ments velocities  and  accelerations  of  the  nodes. 

Equations  (2.1)  to  (2.4)  can  be  written  in  indicial  notation  as 


A.  = A.  .A. 
1 1]  D 

A.  = A. .A. 
1 Di  D 

A.  = y. .A. 
1 ID  D 

A.  = y . . A. 
1 Di  D 

The  coordinate  system  in  which  a vector  is  expressed  will  hence- 
forth be  designated  by  the  bars  and  hats.  Thus  the  components  of  a 
vector  K in  terms  of  the  body  coordinates  of  rigid  body  I are  denoted 
by  A^j,  i = 1 to  3 denoting  A^j f ^yl'  ^zl'  ^^spectively.  Furthermore, 
the  set  of  three  Cartesian  components  is  often  written  as  a matrix 
as  in  Eq.  (2.3). 

The  orientation  of  a node  is  described  by  the  unit  vect:ors 
^il'  '^hile  the  angular  velocities  and  angular  accelerations  are 
treated  in  body  coordinate  components,  and  respectively. 
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The  forces  and  moments  at  the  nodes  are  similarly  denoted  by 

and  respectively,  and  may  be  subdivided  into  externally  applied 

forces  and  moments, F.,  and  M.^  » and  the  forces  and  moments  due  to 

.lx  11 

the  resistance  of  the  deformable  elements,  F^”^and 

il  il 

Only  two  unit  vectors  for  the  nodal  coordinates,  bj^j.  and  b^j, 
are  stored  per  primary  node.  The  third  unit  vector  is  then  found  by 


(2.6) 


This  method  thus  employs  six  numbers  (three  components  of  two  vectors) 
to  describe  the  three  rotational  degrees  of  freedom.  Though  this  at 
first  appears  somewhat  wasteful,  it  should  be  noted  the  alternative, 
a description  by  Euler  angles,  has  serious  shortcomings: 

1.  Euler  angle  formulations  are  not  linearly  independent  for  all 
values  of  the  Euler  angles. 

2.  The  generalized  moments  corresponding  to  Euler  angles  are  not 
easily  interpreted  in  a physical  sense. 

3.  The  equations  of  motion  and  the  transformations  between  body  and 
global  coordinates  in  terms  of  Euler  angles  are  complex  and  computa- 
tionally demanding  because  they  involve  many  trigonometric  functions. 

All  six  components  are  stored,  because  if  only  a total  of 
three  of  the  six  components  of  the  two  unit  vectors  were  stored 
(with  the  remaining  components  computed  from  the  fact  that  the 
two  body  vectors  are  orthogonal  and  unit  vectors) , then  the  re- 
maining components  would  have  to  be  determined  from  square  roots. 

The  signs  of  these  components  could  not  readily  be  determined. 

Element  nodal  forces,  moments,  displacements  and  rotations  are 
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denoted  by  f^j»  "'il‘  *^ii'  ®il'  ^®spectively , where  I denotes  the 
generic  node  number,  which  for  each  element  ranges  from  1 to  the  num- 
ber of  nodes  in  the  element.  Sometimes  a superscript  is  used  to 

• 1 • • . G ^ 

indicate  the  pertinent  element,  i.e.  force  components  at 

node  I of  element  e. 

The  inertial  properties  are  described  by  the  translational  masses 
of  the  primary  nodes,  p^,  and  the  principal  moments  of  inertia  of 
the  primary  nodes,  ^yyj'  ^zzJ*  angular  momenta  of  the 

nodes  are  then  given  by 


L . _ * I . , ,0),  _ 
3J  3kJ  kJ 


(no  sum  on  J) 


(2.7) 


The  element  quantities  are  extracted  from  the  global  quantities 
in  the  usual  manner  by  a Boolean  connectivity  matrix  l'  , so  that 

A1 


(e)  , (e)  „ 

“iA  = “il 


(2.8) 


where 


=1  if  the  Ath  generic  node  of  element  e corresponds 
to  the  Ith  primary  node  of  the  system. 


= 0 otherwise 
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2.  Derivation  of  Equations  of  Motion 

We  consider  here  the  development  of  the  equations  of  motion  for 
the  assemblage  of  rigid  bodies  and  deformable  elements.  The 
equations  are  obtained  from  the  principle  of  virtual  work  with  the 
inertial  forces  included  in  a d'Alembert  sense.  The  principle  of 
virtual  work,  when  applied  to  the  system  treated  here,  states  that 


Me)*-(e)  .•'(e)*  ;;(e) 

'"iA  ^iA  “iA  "^iA 


il  il  il  il 


(2.9) 


- Pi  ^il*  ^il  - “il^il 


where  the  superscript  e is  summed  over  all  elements,,  subscript  A over 
all  nodes  of  each  element,  and  I over  all  primary  nodes.  Superscript 
dots  denote  time  derivatives,  while  asterisks  denote  virtual  quantities 

The  left  hand  side  of  Eq.  (2.9)  represents  the  rate  of  work 
expended  on  the  deformable  elements,  that  is,  the  internal  rate  of 
work,  while  the  first  two  terms  of  the  right  hand  side  represent  the 
rate  of  external  work.  The  rate  of  work  of  the  inertial  forces  is 
represented  by  the  last  two  terms  of  the  right  hand  side. 

To  obtain  the  equations  of  motion,  the  virtual  nodal  velocities 
of  the  element  on  the  left  hand  side  of  Eq.  (2.9)  must  be  expressed  in 
terms  of  the  global  virtual  nodal  velocities.  In  deriving  these 
expressions  we  will  separately  consider  the  case  when  A is  a primary 
node  and  when  A is  a secondary  node. 
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When  A is  a primary  node,  the  required  relationships  are  obtained 
directly  from  Eq.  (2.5),  which  gives 


y . . . 

]i  AI  I] 


(2.10) 


le)  - 

u . . A i.'  0),  _ 

]k  AI  kl 


(2.11) 


ii  I 


j 


When  A is  a secondary  node,  the  required  relationships  are 
developed  as  follows.  We  note  that  whenever  A is  a secondary  node 
associated  with  a primary  node  J,  then  both  nodes  are  points  of  a 
single  rigid  body,  so  that 


- (e) 

w . ' = U) . _ 

lA  ? I 


(2.12) 


and  consequently  Eq.  (2.11)  follows  for  the  angular  velocities.  To 
obtain  the  counterpart  of  Eq.  (2.10),  we  first  designate  the  vector 
from  I to  A by  Because  both  points  are  on  the  same  rigid  body, 

the  components  of  this  vector  in  the  body  coordinates  will  not  vary  with 
tiri.c . The  global  components  of  this  vector  are  given  by 


= A . .X_, 
iIA  1]  ]IA 


(2.13) 


while  the  global  components  of  the  initial  vector  between  I and  A 
are  given  by 


O O . 

x.^-  = A . . X.__ 
iIA  1]  ]IA 


(2.14) 
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where  the  superscript  nought  denotes  the  original  (i.e.  at  time  zero) 
value  of  the  variable.  The  displacement  of  the  secondary  node  is 
then  given  by 


’'‘iA 


X.*  - 
lA  xA 


X._  + X._- 
xl  iIA 


(2.15) 


^il 


^ilA 


By  substituting  Eqs.  (2.13)  and  (2.14)  into  Eq.  (2.15),  we  find 


'^iA  ” '^il  ^ ^ij^ilA  ^ij^jlA 


(2.16) 


Transforming  the  above  to  the  body  coordinates  of  body  I and  talcing 
its  derivative  with  respect  to  time,  we  find 


u.,  = u._  + il..  aj._ 
xA  Xl  X] 


(2.17) 


where 


n. . = 

ID 


0 

^lA 

"^lA 

lA 

0 

*'lA 

lA 

“’^lA 

0 

(2.18) 


By  again  applying  the  appropriate  transformations  from  Eqs.  (2.5), 
we  then  find 


'^iA^  “ ^ji*'AI^"jI  ^ji^jAjl*'AI^“£I 


(2.19) 
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equa- 


Egs.  (2.10),  (2.11),  and  (2.19)  hold  both  for  the  actual 
velocities  and  the  virtual  velocities,  if  we  substitute  these 
tions  into  the  left  hand  side  (LHS)  of  Eq.  (2.9),  using  Eq.  (2.10) 
whenever  node  A is  a primary  node  and  Eo.  (2.19)  whenever  A is  a 
secondary  node,  we  obtain 


LHS  = u*jPjjt  ^ -*^j^nt 


(2.20) 


where 


and 


pint 


AI  3A 


5 (e)- (e) 
AI  kA 


Ad  ^ji  iA 


(2.21) 

(2.22) 


(2.23) 


for  both  primary  and  secondary  nodes  A;  while  if 


A is  primary 


(2.24) 


and  if  A is  secondary 


(2.25) 


Equations  (2.23)  to  (2.25)  may  be  written  in  matrix  form  as  follows 

{f^®b=  []i]{h 


(2.26) 
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} = [ X]'^ [y]  {m} 


(2.27) 


= [X]’’[y]{m}  + [Ql’’lX]’’[y  ] {f } 


(2.28) 


Thus  for  primary  nodes,  the  nodal  forces  and  moments  are  simply 
related  by  the  coordinate  transformations,  while  for  secondary  nodes 
an  additional  moment  is  introduced  in  the  transformation  because  of 
the  moment  arm  effected  by  the  vector  between  the  nodes.  The  total 
internal  nodal  forces  are  obtained  from  the  terms  given  in  Eqs.  (2.23) 
to  (2.25)  by  Eqs.  (2.21)  and  (2.22).  The  latter  equations  just  re- 
present an  appropriate  summation  of  the  element  forces,  for  as  can 

(e) 

be  seen  from  Eq.  (2.8),  are  Boolean  matrices  consisting  of  ones 

and  zeroes. 

The  equations  of  motion  are  now  obtained  by  substituting  Eq. 
(2.20)  into  Eq.  (2.9),  which  gives  (after  a change  of  dummy  indices 
and  collection  of  terms) 


a *(F^f 

xl  xl 


_ext  , ..  . 

Fii  t PjU.j) 


. ~ * /Hint 
+ 03.  _ (M.  _ 

xl  xl 


(2.29) 


Since  the  virtual  velocities  are  arbitrary,  the  terms  within  the 
parentheses  must  vanish.  The  terms  within  the  first  parenthesis 
immediately  yield  the  translational  equations  of  iriotion 


pa.  = 

Pi  il  ‘'il 


(no  sum  on  I) 


(2.30) 


.A 


We  note  from  Eq.  (2.7)  that 

^jJ  " ^jkj“kJ  ®ij£^ikj“kj“£j  J)  (2.31a) 

where  e. is  the  alternative  tensor.  By  noting  that  the  quantity 

in  the  second  parenthesis  must  vanish  because  of  the  arbitrariness  of 
tojj. , we  obtain 


jJ  "jJ 

(no  sum  on  J) 


(2.31b) 


Since  the  coordinates  are  prinicpal  coordinates  of  J..,  we  can 
write  these  equations  as  follows 


+ lljzi  - iyyl'“yi“zl  ' "xf  ‘ »xf 


'yyi“yi  * <‘xxi  ' 'zzi'“xi“zi  = - fi- 


= «ext  _ -int 


I ^ + (I  - T )/Ti  M - nint 

ZZI  Zl  yyl  -^xxl^“xl“yl  “ ”zl  " ”zl 


(2.32) 


(no  sum  on  I) 


These  are  the  rotational  equations  of  motion,  and  they  correspond 
to  the  Euler  equations. 


3.  Deformable  Elements 


The  deformable  elements  are  treated  by  a rigid-convected  (or 
corotational)  formulation  previously  described  by  Belytschko  and 
Hsieh  (1972) . In  this  technique,  the  displacements  of  each  element 
are  decomposed  into  rigid  body  displacements  r^  and  deformation 
displacements  d^ 


u . 

1 


r . 

1 


+ d. 
1 


The  strains  are  then  given  by 


2\3xj  3xi/ 


(2.3 


(2.3 


Belytschko  and  Hsieh  have  shown  that  the  matrix  cor- 
responds to  the  difference  between  the  stretch  tensor  and  the 
unit  tensor.  Hence,  this  strain  corresponds  closely  to  the  common 
definition  of  engineering  strain. 

For  purposes  of  developing  element  relations,  Eq.  (2.33)  is 
expressed  in  matrix  form 


{e}  = (E]{d}  (2.3 

where  {d}  is  the  matrix  of  nodal  deformation  displacements.  It  is 
also  necessary  to  find  a matrix  (T]  such  that 


It  then  follows  that  the  nodal  element  forces 

{f^®h  = tT]'^{f‘^} 

{f*^}  = f[E]'^{a)  dV 
V 


} are  given  by 


(2.37) 

(2.38) 


where  V is  the  volume  of  the  element,  {0}  the  stresses  measured  in 
the  corotational  coordinates,  and  {f*^}  the  nodal  forces  conjugate  to 
{d},  so  that 


= {d}'^{f‘^}  (2.39) 

where  is  the  internal  work.  Note  that  both  {e}  and  {a}  are 

measured  in  corotational  coordinates,  so  their  rates  or  increments 
are  frame  invariant  and  may  be  used  directly  in  incremental  constitu- 
tive equations  without  any  corrections  for  rotations. 


Spring  Element.  Consider  a spring  element  with  nodes  I and  J.  The 
deformation  of  the  spring  is  completely  defined  by  its  change  in  length 

ijj  = (2.40) 

where  I and  1°  are  the  current  and  original  lengths  of  the  element. 
Since  direct  use  of  this  formula  will  result  in  large  round-off  errors, 
an  alternative  formula  was  used  which  is  derived  as  follows.  If 
the  displacement  and  position  of  the  nodes  are  cor  idered  to  be  vectors. 
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{x}j  - {x}.j.  = {x}°  - {x}^  + {u}j  - {u}j 


Taking  the  scalar  product  of  each  side  of  this  equation  with  itself 
if  follows  that 


t 2{x}^^^{u}j,  t {u}j/{u}j. 


where 


{x}ji  = {x}j  - {x}^  etc. 


after  rearranging  and  factoring,  it  follows  that 


i-i°  = — i-  [2{x}  + {u}^/{u}^^] 

u J.  *J  X wX  JX 


or  ir.  component  form 


6__  = l-SL  = [2(x^_u  tt  + z__u 

IJ  .,.0  JI  xJI  JI  yJI  JI  zJI) 


+ u^  +u^  +u^  1 

xJI  yJI  ZJI^ 


The  element  strain  is  then  given  by 


The  nodal  force  conjugate  to  6 is  the  tension  T = f hence  by 


Eq.  (2.38) 


(2.45) 


= AO  (2.46) 

where  A is  the  area  of  the  cross-section. 

.N  A 

The  stress  o may  be  an  arbitrary  function  of  e.  In  this  program, 
the  stress-strain  law  is 

0 = + k2e’  + 26  t (2.47) 

where  k^^,  and  6 are  constants  input  b/  the  user;  k^^  is  the  linear 
spring  constant.  The  last  term  is  a linear  viscosity  with  6 the 
fraction  of  critical  damping  for  the  vibration  of  this  element. 


Beam  Element.  Consider  a generic  beam  element  with  nodes  I and  J us  shown  In 
Fig.  11.  The  A axis  always  moves  with  the  beam  element  so  that  it 

A 

]oins  the  two  nodes,  while  the  y axis  is  considered  to  rotate  with 
the  beam  in  the  sense  that  its  rotation  is  an  average  of  the  rotation 

A 

of  the  two  nodes  about  the  x axis.  The  deformation  of  the  beam  is 
then  defined  by  its  displacements  relative  to  the  rigid-convectad 
coordinates  (x,y,z)  of  the  element.  The  deformation  displacements 
are  given  by 
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V'  V' 


(2.48) 


where  6j.j  = elongation,  computed  by  Eq.  (2.44) 

/\ 

= torsional  deformation 

^ y«.  A A 

= bending  deformation  rotations 

yl  yJ  zl  zj 

Because  of  the  way  the  motion  of  the  x-axis  is  defined,  the 
rigid  body  motion  r^  need  not  be  computed  explicitly  for  the  de- 
composition needed  in  Eqs . (2.33)  and  (2.34):  the  quantities  de- 
fined above  define  the  deformation  of  the  element  directly  regard- 
less of  the  extent  of  rigid  body  rotation. 

For  the  purpose  of  computing  the  relative  rotations,  0^^.,  0^^, 
0yj,  0^j,  and  the  body  components  of  the  unit  vectors  e^“  and 

02“  (superscript  noughts  denote  the  vectors  in  the  undeformed  con- 
figuration) must  be  stored  for  each  of  the  two  nodes  of  the  element. 
From  a knowledge  of  the  body  components  of  and  , the  rotations 
are  found  as  follows.  Since  the  vector  e^^”  rotates  with  the  node,  it 
indicates  the  direction  of  the  axis  of  the  ele  lent  if  there  was  no 
deformation?  the  angle  between  and  e^^  indicates  the  magnitude  of 
the  deformation,  this  is  illustrated  in  Fig.  11.  Thus 


^ ^ 

6 e„  + 0 e.  = e,  x e,' 
y 2 z 3 1 1 


(2.49) 


For  the  purpo:  j of  this  computation,  we  transform  the  components  of 
e^®  from  body  to  element  coordinates  using  Eqs.  (2.3)  and  (2.4)  so 


1 e® 

1 ®lx 

1 ^ 

. , T 

e® 

lx 

< e? 

1 

= [p]  [X] 

e? 

ly 

U® 

'®lz 

e° 

®lz 

Then  substituting  into  Eq.  (2.49),  we  find 


« ^ ->■ 

0 e_  + 6 e. 
y 2 z 3 


= det 


'1 
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3 

0 


L'=lx 


'ly 


IzJ 


-e,  e„  + e,  e. 
Iz  2 ly  3 


(2.51) 


Thus 


6 = 
y 


'Iz 


0 = e? 

z ly 


(2.52) 


The  deformation  torsional  rotation  is  found  by  taking  the  cross- 
product  of  e°j  and  e®j  and  projecting  this  vector  on  the  current 
axis  of  the  beam.  This  yields 


xIJ 


^®2i  ^ ®2j^  = 


-f 

->• 

®2 

-*■  - 
e 

3 

0 

x2I 

e° 

y2i 

"z2I 

x2J 

e° 

y2J 

z2J_ 

'y2I  z2J  y2J  z2I 


(2.53) 


Eqs.  (2.52)  and  (2.53)  require  the  assumption  that  the  deformation 
displacements  of  an  element  be  small.  This  implies  that  the  relative 
rotations  are  sufficiently  small  so  that  the  decomposition  of  the 


rotation  vector  into  0^  and  0^  implicit  in  Eq.  (2.49)  be  valid.  How- 


ever, the  overall  rotation  of  the  beam  element  may  be  arbitrarily 
large . 

The  deformation  displacement  field  for  the  beam  element  is  con- 
sidered to  consist  of  transverse  displacements  that  are  cubic 


functions  of  x,  while  the  axial  displacement  is  a linear  function  of 


X.  This  can  be  written 


dx  = (1-C) 


(2.54a) 


dy  = + (-52+53)^0 


(2.54b) 


(2.54c) 


where 


6 = £9 

X ^ xJI 


(2.54d) 


£=  ^ 
^ i 


(2.54e) 


and  X is  taken  to  originate  at  node  I;  the  superscript  m is  used 
to  indicate  that  these  are  the  displacements  of  the  mid  surface. 

If  we  impose  the  usual  Euler-Bernculli  beam  assumptions  that  normals 
to  the  midline  remain  straight  and  normal,  we  obtain 


a = d"*  - y — Z 

V V J 


- . '‘“x 


(2.55a) 


(2.55b) 


d^  = d-  + ye. 


(2.55c) 


where  H(y,z)  is  the  warping  function. 

The  strain  displacc-ment  equations,  Eq.  (2.34),  can  now  be  written 
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These  equations  are  valid  as  long  as  {3d^/3x)*  and  (3d^/3x)*  are 
small  compared  to  (3d^/3x)^.  Although  this  condition  is  similar  in 
appearance  to  that  of  moderate  rotation  theories,  it  is  far  less 

A A 

restrictive,  because  d^  and  d^  are  the  displacements  relative  to  the 

A A 

corotational  coordinate  x.  By  reducing  the  size  of  the  element,  d 


and  d^  can  be  made  as  small  as  reeded,  albeit  at  a rapidly  increasing 


cost. 


For  the  beam 


and 


A T'  A A A 

{e}  = (e  , 2c  , 2e  ) 
' X xy  xz 


[E]  = 


z(6C-4)  y(4-65)  z(6C-2)  y(2-6C) 


0 (3H/3y)-z 

0 (3H/3z)-y 


0 

0 


0 

0 


0 

0 


(2.57) 


(2.58) 


I w 

'I 


So,  having  defined  the  strain-displacement  relations,  the  stress 

{o}^={a,o  ,cr  } can  be  computed  for  any  constitutive  law,  and 
X xz 

then  the  deformation  nodal  forces  by  Eq.  (2.38). 

Referring  to  Eqs.  (2.39)  and  (2.40),  it  may  be  noted  that  {f^} 
for  the  beam  is  given  by 
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(2.59) 


(f  ^iiTi  •■in  _fin  _/in  _) 

xJ  xJ  yl  zl  yJ  zJ 


The  complete  nodal  force  matrix  could  be  computed  from  Eq.  (2.37) 
after  [T] , as  defined  by  Eq.  (2.36),  is  found.  However,  it  is 
simpler  to  first  find  {f}  by  equilibrium  and  then  transform  to  '>e 
global  components  by  a simple  vector  transformation;  the  two  pro- 
cedures are  equivalent.  Equilibrium  yields 


f T = -f  T 

xl  xJ 


m _ = -m  , 
xl  xJ 


'zJ 


m _ + m , 
yi YjL 


^zl  " "^zJ 


"'zl  + ""zJ 


yJ 


^yl  - ‘^yJ 


(2.60) 


The  global  components  can  then  be  found  by  Eq.  (2.1)  and  (2.2). 

As  indicated  in  Chapter  I,  the  program  has  two  methods  for 
treating  the  moment  stiffness  of  the  beam.  For  linear  beams  or  beams 
with  simple  nonlinear  characteristics  where  the  overall  bending 
moments  can  be  determined  as  functions  of  the  deformation  rotations, 
the  bending  moments  may  be  determined  directly  by  Eqs.  (1.4). 

The  second  method  for  obtaining  the  nodal  moments  is  by  inte- 
grating the  stresses  appropriately  through  Eqs.  (2.38),  with  [E]  de- 
fined by  Eq.  (2.58).  This  integration  is  done  only  to  obtain  the 
axial  force  and  moment.  The  numerical  integration  for  torsion  is 
included  neither  in  this  description  nor  in  the  program.  If  torsion 
were  included,  both  the  state  of  stress  and  strain  would  be  biaxial, 
which  increases  computations  immensely.  Since  Euler-Bernoulli  theory 
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does  not  include  the  shear  stresses  that  arise  from  bending,  only 
the  shear  strains  due  to  torsion  are  available.  Thus  the  additional 


computations  would  not  appear  to  be  justified. 

The  integral  for  the  nodal  forces,  Eq.  (2.38),  is  evaluated 
numerically.  Many  investigators  have  used  Gaussian  quadrature 
formulas  for  this  integration.  However,  these  formulas  involve  far 
more  computations  than  trapezoidal  integrations,  with  little 
advantage  in  accuracy.  The  only  drawback  of  a strict  trapezoidal 
integration  is  that  the  elastic  bending  stiffness  of  the  element  is 
not  recovered.  A modified  trapezoidal  integration  originally 
described  in  Marchertas  and  Belytschko  (1974)  is  used  here.  The 
beam  is  subdivided  into  strips  as  shown  in  Fig.  3.  In  each  strip, 
as  for  example  the  strip  indicated  by  cross-hatching  in  Fig.  3,  the 

A 

axial  stress  o is  assumed  to  vary  bilinearly  in  n and  so 


o(n,5)  = (1-5)  (1- n)  + 0j^n(i-5) 

" -.j(l-r,)5  . G.j  r,  5 


(2:6l)i 


Substituting  this  form  into  Eq. 
5,  we  find  the  expressions  for 


(2.37)  and  integrating  over  r)  and 

''  j ^ P 

f*''-  and  m , to  be 


^ = 7 th  [a._  + o._  + 0.,  + o.  ] 
xJ  4 il  3 1 iJ  3J 


(2.62) 


m^^  = i th  [2o.  z.  + 2o._z.  + o.  z . + o.^z.] 
yJ  6 ij  1 ]J  3 ij  3 3J  1 


(2.63) 


The  formulas  for  m m , and  m ^ are  similar.  The  total  nodal  forces 

yl  zJ  zl 


are  then  found  by  adding  the  contributions  of  all  of  the  layers.  The 


major  distinction  from  a trapezoidal  integration  is  in  Eq.  (2.63), 

/V 

which  is  based  on  a form  quadratic  :.n  z.  This  modification  preserves 
the  exact  elastic  bending  stiffness  of  the  element.  Moreover,  the 

A 

assumption  of  a linear  stress  variation  in  x implies  a linear  moment 
field,  which  is  consistent  with  equilibrium  of  the  element.  The 
approximation  errors  therefore  lie  entirely  in  the  stress-strain  law, 
which  is  necessarily  satisfied  only  at  the  end  points  of  the  element 
for  nonlinear  materials. 


Hydrodynamic  Element.  The  hydrodynamic  element  is  a pentahedron  with 
six  nodes  at  its  corners,  as  shown  in  Fig.  4.  The  triangular  sur- 
faces of  the  pentahedrons  are  called  the  top  and  bottom  surfaces, 
respectively,  and  each  of  these  is  defined  by  three  nodes.  All  of 
the  nodes  on  each  of  these  surfaces  must  be  associated  with  the  same 
primary  node,  for  the  top  and  bottom  surfaces  are  considered  rigid. 
The  generic  element  is  nuiobered  counterclockwise  and  viewed  from  top 
to  bottom.  The  edges  of  the  element  arc  defined  by  the  lines  connect- 
ing nodes  I - L,  nodes  J - M,  and  nodes  K - N,  respectively.  The 
deformation  of  the  hydrodynamic  element  is  characterized  by  the 
dilatation.  A,  which  is  given  by 


V - V' 


where  V is  the  volume  of  the  element,  and  a superscript  nought  denotes 
the  oriainal  volume.  For  the  purpose  of  computing  the  volume  and 
change  in  volume,  the  pentahedron  is  subdivided  into  three  tetrahedrons . 
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The  volume  of  any  tetrahedron  vrith  generic  node  numbers  I,J#K  and  L 
is  given  by 


^JI 

^KJ 

^LK 

''lOKL  - ? 

yji 

^KJ 

^LK 

2 

z 

2 

L Ji 

KJ 

LK 

6 ^^JI^^LK^KJ  " ^KJ^LK^  ^KJ^^JI^LK  “ ^KL^JI^ 


^IJ  = “ ^J' 


(2.65) 


The  volume  of  the  pentahedron  is 


V = V + V + V 

IJKL  JLMN  JKLN 


(2.66) 


To  define  the  matrix  [E]  of  Eg.  (2.35),  we  ta)ce  the  time  derivative 
of  the  dilatation  as  defined  in  terms  of  the  nodal  coordinates  in 
Eqs.  (2.66).  Noting  that  the  time  derivatives  of  the  nodal  coordin- 
ates are  the  velocities,  we  obtain  [E] , which  is  given  in  Table  1. 

The  pressure  in  the  hydrodynamic  element  is  assumed  to  be  con- 
stant, so  that  the  integral  expression  for  the  nodal  forces  given  by 
Eq.  (2.38)  yields  that  the  nodal  forces  are  giver,  by 

T 

{f}  = -P  (El  (2.67) 
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^ 'fc* . ..  * *>'*i>*"  T ijv'':  V* 


■3 

■J 


Table  1.  The  [E]  Matrix  for  Hydrodynamic  Elements 
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The  pressure  is  computed  by 


p = BA  + aBA  (2.68 

where  a is  a linear  viscosity  coefficient  and  B is  a bulk  modulus. 

For  the  nodes  that  are  secondary,  the  nodal  forces  are  then  trans- 
formed by  Eq.  (2.28). 


Elastic  Planes.  The  ejection  seat  geometries  are  represented 
by  an  assemblage  of  piecewise  linear  surfaces,  which  may  be  either 
elastic  or  rigid.  This  assemblage  of  planes  may  translate  in  an 
arbitrary  direction  in  space  as  an  arbitrary  function  of  time,  which 
is  prescribed  in  the  input.  When  any  of  the  rigid  bodies  contacts 
one  of  the  planes,  its  equations  of  motion  are  modified  so  as  to 
reflect  its  interaction  with  the  plane. 

Each  of  the  piecewise  linear  surfaces  is  defined  by  three  points 
X®,  I - 1 to  3,  in  the  original  configuration.  The  coordinates  of 
these  points  are  given  in  the  global  system.  The  three  points  must 
be  oriented  so  that  the  normal  vector  as  obtained  by  the  right  hand 
rule  points  in  the  direction  of  the  skeletal  model.  In  order  to 
establish  the  required  equations,  the  normal  to  the  surface  it  first 
computed  by  the  formula 


n 


(2.69 


The  minimum  distance  between  any  node,  x , and  the  surface  is  given 


D“ 


p • n 


(2.70) 


where  p®  is  a vector  from  any  point  in  the  surface  to  x . As  the 

P 

surface  translates,  (rotations  of  the  surface  are  not  included) , the 
unit  normal  remains  unchanged,  while  the  positions  of  and  x^ 
become 


(2.71) 

x = X®  + u 
p p p 

where  u^  is  the  prescribed  displacement  of  the  ejection  seat.  The 
distance  of  point  P from  the  surface  can  then  be  computed  by 

D = p • n (2.72) 

When  the  distance  of  the  node  from  the  surface  is  negative  or  zero, 
the  node  is  assumed  to  be  in  contact  with  the  surface.  It  is  tiien 
necessary  to  check  its  acceleration  normal  to  the  plane.  If  the 
acceleration  is  into  the  plane,  the  node  will  remain  in  contact  with 
the  plane,  whereas  when  the  acceleration  is  away  from  the  plane, 
there  is  no  contact  and  no  modifications  of  the  equations  of  the 
node  are  needed. 

If  the  node  is  in  contact  with  the  plane,  a force  normal  to 
the  plane  is  applied  to  the  node  v’hich  is  proportional  to  the 
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CHAPTER  III 


MATHEMATICAL  ASPECTS  OF  SOLUTION  PROCEDURES 


1.  Explicit  Time  Integration 


One  of  the  i.ethods  for  solving  the  equations  of  motion  of  this 
model  is  the  explicit  method.  This  method  is  by  far  the  most  effici- 
ent per  time  step,  but  the  time  step  must  be  quite  small  if  the  model 
has  a high  frequency  content,  that  is,  low  masses  connected  by  stiff 
springs.  The  stability  limits  are  discussed  subsequently  in  this 
Section. 

The  explicit  integration  employs  the  Newmark  S“formulas  (Newmark 
(1959))  with  6=0,  which  are  almost  identical  to  the  central  difference 
formulas  (see  Belytschko  (1974)).  These  formulas  predict  the 
velocities  and  displacements  at  a time  step  in  terms  of  the  accelera- 
tions at  the  previous  step.  For  the  translational  components,  these 
formulas  may  be  used  directly,  so 


u?""^  = u?  + i At(u?  + u?-^^) 
il  il  2 ' il  il  ' 


(3.1) 


“n' = “II  * “ '^ii  * J “il 


(3.2) 


where  the  superscripts  denote  the  time  step  and  At  the  time  increment 
during  a step. 

For  the  rotational  degrees  of  freedom,  Eqs.  (3.2)  cannot  be  used 
directly  because  the  orientations  are  described  by  unit  vectors 
b^^  and  its  rates  are  not  equivalent  to  the  angular  velocities  and 
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accelerations.  The  counterpart  of  Eq.  (3.2)  is 


db?  -I  ^2^3 

+ A t i At 

^ dt  2 ^^.2 


2 1 


From  vector  analysis 


db . 

— i = S X s. 

dt  ^ 


(3.3) 


(3.4) 


i-  = 0)  X (o)Xb.)  + (axb.)  (3.5) 

dt^  ^ "■ 

and  substituting  Eqs.  (3.4)  and  (3.5)  into  Eq.  (3.3)  yields  a formula 
for  the  updated  unit  vectors 

-►i+1  ->“i  ->-i 

b.  = b.  + At  (loxb . ) 

1 1 1 

+ ^ At^  [io  X ((oxb|)  + (axb|)  3 (3.6) 


To  obtain  a formula  for  the  updated  components  of  the  unit  vectors, 
we  temporarily  fix  the  x^  coordinates  and  consider  a particular  unit 
vector  b^  in  Eq.  (3.6);  we  then  dot  Eq.  (3.6)  with  the  unit  vector 
corresponding  to  the  component  of  b^  to  be  updated.  For  example,  the 
updated  x-component  of  is  found  by  letting  i-3  in  Eq.  (3.6)  and 
talcing  the  scalar  product  of  both  sides  of  the  equation  with  b^^, 
which  yields,  after  some,  simplification 
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(3.7a) 


A (w^oi^+a^  ) 
2 X z y 


Similarly 


' b^^^  = ^ to)^  + i A t^  (w^w^-a^)  (3.7b) 

jy  /,  -j  X ^ y z X 

b|y^  = b^  * b^"^^  = Ato)^  + I A tMoj^oJ^+a^)  (3.7c) 


Normality  and  orthogonality  are  used  to  find  the  remaining  com- 
ponents of  and  and  it  is  assumed  that  the  rotations 

during  a time  step  are  small  so  that  second  order  terms  can  be  neg- 
lected. From  the  normality  of  it  follows  that 


3z 


= [- 


(3.8) 


while,  if  it  is  first  assumed  that  ~ 1,  orthogonality  yields 


3x  ly  3y  3z 


(3.9) 


The  component  can  then  be  found  by  normality 

The  components  of  the  b^^  and  b^  vectors  given  above  are  in  terms 
of  the  coordinates  of  that  node  at  time  step  j.  Using  Eq.  (2.3) 
with  the  [A]  matrix  defined  in  the  previous  time  step,  the  components 
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are  transformed  to  global  coordinates;  the  vector  is  then 

found  by  a vector  product. 

The  tlowchart  for  this  numerical  procedure  is  shown  in  Table  2. 
Besides  the  approximations  of  small  rotations  implicit  in  Egs.  (3.9) 
and  (3.10),  the  major  source  of  truncation  error  lies  in  the  evalua- 
tion of  the  angular  accelerations  from  the  equations  of  angular 
motion,  (2.32).  These  equations  involve  quadratic  combinations  of 
the  angular  velocities  at  the  same  time  step  as  angular  accelerations, 
so  that  the  new  angular  velocities  must  be  approximated.  This  source 
of  error  can  be  mitigated  by  an  iterative  procedure,  i.e.  finding 
the  angular  accelerations  from  an  approximation  to  Eq.  (2.32),  inte- 
grating these  accelerations  to  find  the  new  angular  velocities,  and 
then  repeating  the  computation  of  the  angular  accelerations.  This 
procedure  was  tried,  but  it  was  found  to  have  insignificant  effects 
on  the  solution. 

Another  consequence  of  the  existence  of  the  quadratic  combina- 
tion of  velocities  in  Eq.  (2.32)  is  that  the  classical  Fourier  theories 
of  numerical  stability  are  not  applicable  even  for  linear  material 
problems.  Therefore,  it  has  been  necessary  to  pick  time  steps  with 
the  standard  formula  in  terms  of  the  maximum  eigenvalue  X of  the 
system 


At^  <f 


only  serving  as  a guideline.  To  insure  post  facto  that  the  computa- 
tions are  stable,  energy  balance  checks  were  made  as  follows. 
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Table  2 

Flow  Chart  of  Explicit  Integration  Procedure 

Set  initial  conditions,  t=0 
Compute  {u(t+At)}  by  Eq.  (3.2) 

Update  unit  vectors  bj  by  Eqs.  (3.7)  to  (3.10)  and  transform 
to  global  components  by  Eq.  (2.3) 

Find  deformation  displacements  {d}  by  Eqs.  (2.44),  (2.49)  and 
(2.53) 

A 

Find  the  strain  in  the  convected  coordinates,  (e)  by  Eq.  (2.34) 
Stress-strain  law 

Find  local  nodal  forces  {f^}  by  Eq.  (2.38) 

Add  {f^}  into 

Compute  {vi(t+At)}  by  Eqs.  (2.30)  and  (2.31) 

Compute  {0(t+At)}  by  Eq.  (3.1) 
t = t + At;  go  to  2 


The  kinetic  energy  at  each  time  step  j is  defined  by 


_j  1 .j2  ^ 1 -T  T -3 

T = ■:r  p.UT^  + — 

2 1 il  2 il  iJ.1  i 


ZJ 

il 


where  1 is  sxmmed  over  the  primary  nodes  in  the  system.  The  external 
work  is  defined  by  a trapezoidal  integration  in  time,  so 


w‘"'°  = 0 


ext.j  ^ ^ 1 j _ „j-l  ext.j  ^ 

2 xl  il  il  il 


The  internal  work  is  similarly  defined  by 


= 0 


^int,j  ^ ^int,j-l  ^ 

2 lA  lA  lA  lA 


Energy  balance  then  requires  that 


yjint,j  + .pj  _ j^ext,j  ^ ^^^^int,j  ^ 


If  e is  greater  than  about  0.05,  a solution  is  considered  unacceptable. 
The  source  of  the  energy  error  may  be  an  arrested  or  incipient 
numerical  instability,  or  excessive  truncation  error. 


68 


MM 


\ I 

! •* 


i 


2.  Implicit  Time  Integration 
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The  second  method  for  integrating  the  equations  of  motion  in 
time  is  the  implicit  method.  This  method  permits  the  use  of  rather 
large  time  steps,  particularly  if  the  high  frequency  aspects  of 
the  response  are  unimportant.  However,  the  method  involves  the  solu- 
tion of  nonlinear  equations  in  each  time  step.  The  nonlinear 
equations  are  solved  by  a Newton-Raphson  procedure,  so  that  a sequence 
of  solutions  to  linear  algebraic  equations  are  used  to  obtain  the 
solution  to  the  nonlinear  equations.  The  linear  equations  are  solved 
by  a Cholesky  decomposition  technique;  hence  the  solution  time  is 
bandwidth  dependent.  Therefore  the  method  becomes  quite  uneconomical 
for  the  complex  models,  where  the  bandwidth  is  usually  quite  large, 
for  in  spite  of  the  possibility  of  taking  large  time  steps,  the 
computational  effort  required  per  time  step  is  so  great  as  to  pre- 
clude its  use.  More  specific  figures  for  solution  times  are  given 
subsequently . 

The  implicit  integration  employs  the  trapezoidal  integratic. 
formulas,  which  correspond  to  the  Newmark  ^-formulas  with  3=1/4 
(Newmark  (1959)).  These  formulas  are 


u^^^  = u^  + At  + jAt^(u^  + 
il  il  il  4 il  il 


(3.16) 


(3.17) 


From  the  above,  it  immediately  follows  that 
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where 


..j+1  4 ,,  . ..  i 4 . j 

u-^_  = (Au  ) - u._  - — 

ll  ^^.2  il  il  II 


.1+1  2 , .1 

u.^  = ■: — Au.^  - ur^ 

il  At  il  il 


Au.  = 

il  il  il 


The  above  equations  must  hold  both  for  increments  in  rotations  and 
translations,  provided  that  the  increments  in  rotations  are  small,  so 

s ‘"ii  - “L 

The  increment  A6  has  no  precise  meaning,  for  it  is  neither  referred 

♦i  -►1+1  . 

to  the  bj  or  the  b^  coordinate  systems,  but  to  some  intermediate 

coordinates.  However  A6  is  not  a basic  variable,  for  orientations 

are  ultimately  described  by  the  body  coordinate  system  so  this 

ambiguity  is  not  important. 

The  equations  for  implicit  integration  are  best  developed  in 
matrix  form,  so  we  define  the  nodal  matrices  {D>,  {D}^  and  {F}^. 
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The  complete  matrix  of  nodal  displacements  is  {D} , which  consists 
of  {D}j,  with  I from  the  number  of  nodes  stacked  vertically  in  order 
in  the  column  matrix. 

The  implicit  solution  procedure  requires  a linearized  estimate 
of  the  internal  forces  at  the  end  of  the  time  step  in  terms  of  the 
incremental  forces.  This  is  provided  by  the  tangential  stiffness 
matrix  ih]  which  gives 


^pj+l,intj  ^ ^pj,rntj  ^ [k]{AD} 


(3.24) 


i } I 


The  tangential  stiffness  here  contains  the  damping  terms;  details 
are  given  in  Appendix  2. 

The  basis  of  the  implicit  method  is  to  solve  the  equations  of 
motion  at  time  step  j+1  directly  in  terms  of  the  displacements.  For 
this  purpose,  the  equations  of  motion  are  expressed  at  j+1  in  terms 
of  the  solution  at  the  previous  step.  The  right  hand  sides  of  the 
equations  of  motion,  Eqs.  (2.30)  and  (2.31)  are  expressed  in  terms  of 
increments  of  displacements  by  means  of  Eqs.  (3.24),  which  gives 


Pj(4u.j  - It 


fp j+1, ext 
4 ■ il 


j+l,int. 
^il  ^ 


(3.25) 


^siKi  - 


+ + At  mpAe,,  + 


^tsSu'ir  %“t  ^ ^ Atm- Ae^  + AB^Ae^) 


_ Mi  /Gj+l»ext  r.j+1/int. 


■ 


We  now  (1)  use  Eg.  (3.24)  to  estimate  the  internal  forces  at  j+1, 

(2)  neglect  the  quadratic  terms  A0^A6^  and  (3)  neglect  any  anti- 
symmetric terms  arising  from  the  rate  of  change  of  angular  momentiam. 
Eg.  (3.25)  and  (3.26)  can  then  be  written 


([Ml]  + -4-  [M2]  + ^ [K]){AD} 


= At[Ml]{P}  + ^ ([M1]{D^} 


+ _ {pj,intj  _ (3^27) 


where 
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(3.28) 
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(3.30) 


= e. 
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Eqs.  (3.27)  are  solved  in  each  time  step.  The  error  is  then  computed 
by 


,err  _ j+l,ext  j+l,int 
il  " ^il  ■ ^il  " 


i^err  ^ -j+l,ext  _ j^j+l,int  _ £j+l 
il  il  il  il 


(3.31) 

(3.32) 


The  error  must  satisfy  the  criterion 

(1D)’^{F^")  < + T^)  (3.33) 

where  e is  a tolerance,  and  and  T^  are  the  internal  and  Icinetic 

energies,  as  defined  by  Eqs.  (3.12)  and  (3.14).  This  criterion  re- 
quires that  the  energy  transfer  to  the  system  arising  from  error  be 
bounded.  If  Eq.  (3.33)  is  not  met,  Eqs.  (3.27)  are  resolved  with 
{F  } added  to  {F  }.  This  does  not  require  much  extra  time,  for 
the  coefficient  matrix  is  triangulated,  at  most,  onlj  once  per  time  step. 

Belytschko  and  SchoeberJe  (1975)  have  shown  that  this 
procedure  leads  to  a stable  algorithm  regardless  of  the  size  of  the 
time  step  as  long  as  Eq.  (3.33)  is  met.  The  cited  proof  was  only  for 
nonlinear  materials,  but  so  far  no  unstable  computations  have  been 
reported,  when  Eq.  (3.33)  was  satisfied  for  a reasonable  tolerance. 

The  unit  vectors  are  updated  as  follows.  The  vector  counterpart 
of  Eqs.  (3.16)  and  (3.17)  are 
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6^*1  . ej  t itej  + ^ (Ej  + Ei*') 


t\*^  = ji  . f (Si  . 5j*i, 


(3.34) 


(3.35) 


Multiplying  Eq.  (3.34)  ^ adding  to  Eq.  (3.35),  we  obtain 


^ - H*" 


and  rearranging  and  using  Eq.  (3.4)  yields 

gj+1  _ f (-3+1 X sj^i)  = sj  + ^ (;s^  X 5j) 


(3.36) 


(3.37) 


To  obtain  a formula  for  the  updated  ccxnponente  of  the  unit 
vectors,  we  temporarily  fix  the  x^  coordinates  and  consider  a parti- 
cular  unit  vector  5^in  Eq.  (3.37).  We  then  take  the  scalar  product 
of  Eq.  (3.37)  with  the  unit  vector  corresponding  to  the  component  of 
to  be  updated.  For  example,  the  updated  x-component  of  b^  is 
found  by  letting  1=3  in  Eq.  (3.37)  and  taking  the  scalar  product  of 
both  sides  of  the  equation  with  which  yields,  after  some  simplifi- 


cation 


»3x^  = 53  • 5^*1  =.  f £J  • (E^^l  X Ei*l)  4 ^ -J  (3.38) 


similarly,  the  y-components  of  b^  is  given  by 


u3+l  Sj  • ^j+1  At  tj  . /-►j+l  :^j+li  At  -j  ,, 

b|y  = b^  b^  “ T ^2  * ^3  ^ ■ T “x  (3*39) 
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The  evaluation  of  triple  scalar  products  in  Eqs.  (3.38)  and  (3.39) 
yields 


?j  • ,‘*‘j+l  ~j+l  v.3+1  “j+1  v.j+1 

bJ  X bj  ) - bJ^  - bJy 


(3.40) 


Hence  Eqs.  (3.38)  and  (3.39)  can  be  written  as 


.j+1  At  /-j+l  .j+1  -j+1  v.j+lx  At  -j 


(2.41) 


3y  2 z 3x  X 3z  ' 2 x 


(3.42) 


Now  if  it  is  assumed  that  1»  Eqs.  (3.41)  and  (3.42)  can  be 

solved  simultaneously  to  yield 


j+1  At  .At  - j+1 j , -j+li  . ,-j  . “j+li 
i..  = -r-  [-S-  tii-’  (ci)-'  + C0-'  ) + (u):'_  + w'’ 


3x 


2 “z  K ^ “x  ’ ^“y  “y 


M J.  At^  -j+1  -j  + lx 
T"  “z  “z  ’ 


(3.43) 


.j+1  At  .At  - j+1 j . -j+lx  /- j . “j+l\i/ 

*’3y  ■ T ‘T  “z  + “y  ) - "J  >1/ 


,,  At®  -j+l~j+l, 
(1  + — :: — 0)  0)  ) 
4 z z 


(3.44) 


In  a similar  manner  the  y-component  of  is  found  to  be 


b^+l  = 

ly 


At 

2 


.At  -j+l/-j  . -j+lv 
1“:^  03  (o)  +0)  ) 

^ 2 X y y ' 


(j3  . 


(1  + 


(3.45) 
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Equations  (3.37),  when  substituted  into  the  above,  yield 


3x 


■At  -j+1  - . \ //I  ^ At^  -j+l-j+1. 

“z  ^ ~r  '^Z  “’z  ^ 


(3.46) 


>3+1  = 

3y 


.At  j-^l,o  ,7  X . At^  -j+l-j+l> 

(-y  U3^  A0y  - A0^)/(1  + — o>J  0)^  ) 


(3.47) 


. j+1  ,At  j+l,o  . .7  X ,,1  . At^  -j+l-j+lx 

‘^ly  = < ^®y  + ~ “x  “x  ^ 


(3.48) 


The  remaining  components  are  then  updated  as  described  in  the  previous 
section. 

Tjie  major  computational  effort  in  this  procedure  is  in  triang- 
ulating Eq.  (3.27).  If  the  number  of  degrees  of  freedom  is  N and 

2 

the  semibandwidth  is  B,  triangulation  requires  NB  multiplications. 

The  semibandwidth  is  given  by 


B = max(I-J)  (3.49) 

where  7 and  J are  the  node  numbers  of  any  primary  nodes  connected  by 
a deformable  element,  either  directlv  or  through  secondary  nodes. 

It  is  thus  of  advantage  to  number  nodes  so  that  B is  minimized. 

The  equations  are  triangulized  at  most  once  in  each  time  step. 
The  iterative  procedure  is  performed  using  the  triangularized  mat- 
rix with  direct  backsubstitution.  This  only  requires  2NB  multipli- 
cations. If  the  system  is  only  midly  nonlinear,  the  equations  are 
solved  with  the  triangular  decomposition  obtained  from  a previous 
time  step.  The  computational  procedure  is  shown  in  Table  3. 
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Table  3 


Flow  Chart  of  Iterative,  Implicit  Integration  Procedure 

1.  Set  initial  conditions:  t=0;  {ii(0)}  must  be  computed  if  initial 
conditions  are  nonzero 

2.  Compute  coefficient  of  {AD}  and  R.H.S.  of  Eq.  (3.27)  with  {R}  = 0 

3.  Solve  Eq.  (3.27)  for  tentative  increment  in  displacement 

4.  Compute  tentative  displacement,  acceleration  and  velocity  by 
Eqs.  (3.21),  (3.19)  and  (3.20),  respectively 

5.  Update  unit  vectors  by  Eqs.  (3.46)  to  (3.48)  and  transform  to 
global  components  by  Eq.  (2.3) 

6.  Find  tentative  strain  by  Eq.  (2.34) 

7.  Stress-strain  law 

8.  Find  tentative  nodal  forces  by  Eq.  (2.38)  and  add  into  internal 
force  array 

9.  Compute  error  force  by  Eqs.  (3.31)  and  (3.32) 

10.  If  energy  error  criteria,  Eq.  (3.33)  in  not  met,  set 

= {F®*^}  + {F®^^>  and  go  to  3 

11.  Solution  for  step  complete;  update  dependent  variables 

12.  t = t + At;  go  to  2 
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3.  Modal  Analysis 


The  modes  of  the  system  are  obtained  by  a linear  eigenvalue 
technique.  The  use  of  this  technique  thus  requires  that  the  material 
properties  be  linearized  in  the  domain  of  interest,  so  that 


= [K]{D} 


(3.50 


The  use  of  a linear  eigenvalue  analysis  also  requires  that  the  pro- 
duct terms  oi  oi  , oi  oi  in  Eqs.  (2.32)  be  neglected:  thus  the  modal 
X z y z 

analysis  requires  that  the  motion  be  restricted  to  two  dimensions 
or  that  the  rigid  bodies  have  equal  moments  of  inertia  about  all  axes 
In  this  study,  the  first  assv;mption  was  usually  made.  In  addition, 
for  purposes  of  simplicity,  the  effect  of  damping  on  the  modal 
behavior  was  neglected.  Provided  that  the  damping  is  proportional 
to  the  stiffness,  which  is  the  case  here,  the  effect  of  damping  on 
the  frequencies  could  then  be  asscertained  after  the  completion  of 
the  analysis  as  described  subsequently. 

The  modal  equations  are  obtained  in  the  usual  manner  (see 
Przemieniecki  (1968))  by  assuming  harmonic  response  of  the  form 


{D(t)>  = {X}e 


(3.51 


Substituting  Eqs.  (3.50)  and  (3.51)  into  (2.30)  and  (2.32)  and  in- 
voking the  above  assumptions,  yields 


( [K]  - 0)^  [Ml]  ){X}  = 0 


(3.52 


which  is  linear  eigenvalue  problems,  with  the  eigenvalues  (u^  cor- 
responding to  the  natural  frequencies,  and  the  eigenvectors  {X} 
corresponding  to  the  natural  modes  of  vibration. 

The  eigenvalue  problem  of  Eg.  (3.52)  is  put  into  standard  eigen- 
value form  as  follows.  Since  the  products  of  angular  velocities 

appearing  in  the  rotational  equations  of  motion  are  not  treated, 

-1/2 

[Ml]  is  diagonal  and  a matrix  [Ml]  ' can  be  constructed  by  taking 
the  reciprocal  of  the  square  root  of  each  diagonal  term  of  the 
matrix.  Defining 


{X}  - [M1]^/^{X} 


(3.53) 


we  obtain  from  (3.52)  and  (3.53)  that 


([M]'^/^[K]  [mT^/^-wMi]){X}  = 0 


(3.54) 


Ilf 


r.i  ! # 


The  eigenvalues  and  eigenvectors  are  found  by  either  of  two  methods : 

i.  a Jacobian  iterat;  ,n  proce’ure  as  implemented  in  subroutine 
EIGEN  in  the  IBM  Scientific  Subroutine  Package  (1967) ? 

ii.  a power  method  iteration  procedure. 

The  first  method  is  relatively  time  consuir. _ng  and  cannot  take 
advantage  of  the  bandedness  of  the  stiffness  matrix.  Hence  it  can 
be  used  only  for  relatively  small  models.  However,  it  provides  all 
frequencies  and  modes  of  the  system. 

The  second  method  can  take  advantage  of  bandedress  and  hence  is 
applicable  to  larger  models.  Howeve.  , in  the  present  version  only  th.. 
lowest  mode  can  be  found.  The  determrnation  of  subsequent  modes 


Vf 


..Vi-' : ■ • •'  ■'.  .♦.'  v«,-ir»‘'£;  ,V. 


would  entail  the  implementation  of  shift  and  sweeping  procedures 
into  the  algorithm;  see  Gourlay  and  Watson  (1973)  or  Belytschko 
(1974)  . 

If  the  system  is  damped  by  a stiffness  proportional  damping  of 
magnitude  a as  used  herein,  then  the  frequencies  of  the  damped  system 
may  be  determined  directly  from  the  undamped  modes.  Standard  har- 
monic analysis  then  shows  (see  for  example,  Hurty  and  Rubinstein 
(1974)),  that  the  frequency  of  the  damped  vibration  is  given  by 

damp  Vi 

w = to  fl  - — ^ — (3.55 
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CHAPTER  IV 


MATERIAL  PROPERTIES  AND  GEOMETRY 

This  chapter  describes  the  geometrical  and  material  property 
data  used  in  the  models  and  the  sources  and  basis  for  this  data.  For 
purposes  of  convenience,  the  description  of  this  data  is  subdivided 
into  sections  dealing  with  the  spine,  the  head  and  pelvis,  and  the 
rib  cage.  The  inertial  properties  are  described  in  a separate  Section. 
Finally,  the  combinations  of  these  data  used  in  models  of  various 
complexity  are  summarized  by  describing  Models  I,  II,  and  III. 

1.  Thoracolumbar  and  Cervical  Spine 

Each  vertebra  of  the  thoracolumbar  spine  is  modelled  as  a rigid 
body.  The  geometry  of  each  vertebra  is  described  by  the  locations 
of  the  secondary  nodes,  v/hich  serve  to  connect  the  various  deformable 
elements  to  the  ertebra.  An  inertial  segment  usually  does  not 
coincide  with  a vertebral  body.  Instead,  each  vertebra  is  embedded 
in  an  inertial  segment,  and  the  primary  node  of  the  segment  is  its 
centroid.  Therefore,  it  is  worthwhile  in  describing  the  geometry  of 
the  vertebrae  to  define  a base  point  which  is  independent  of  the  loca- 
tion of  the  primary  node  of  the  inertial  segment.  This  base  point  is 
chosen  to  be  the  center  of  the  lower  end  plate  of  the  vertebral 
body.  The  geometry  of  each  vertebra  is  then  characterized  by  the 
positions  in  rhe  vertebra's  body  coordinates  of  twelve  additional 
points:  (1)  the  spinous  process  tip;  (2)  left  and  right  transverse 
process  tips;  (3)  the  left  and  right/inferior-superior  articular  facet 
points;  and  (4)  the  left  and  right/inferior-superior  ligamenta  flava 
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points.  A set  of  coordinates  for  vertebra  L3  is  given  in  Table  4. 

This  data  was  obtained  from  a skeletal  model  by  measurement  and  was  ■: 

I 

then  scaled  to  conform  to  the  measurement  of  the  measured  averages 
of  Lanier  (1939)  by  Schultz  and  Galante  (1970) . 

The  second  aspect  of  the  geometry  is  the  positioning  and  orienta- 
tion of  the  vertebrae  in  space.  In  all  of  these  studies,  it  was 
assumed  that  there  is  no  curvature  in  the  frontal  plane,  so  only  the 
curvature  in  the  sagittal  plane  was  included.  The  first  sets  of 
data  were  generated  by  taking  the  previous  model  of  Schultz,  et  al. 

(1974),  which  was  for  a standing  position,  and  reducing  the  curvature 

so  as  to  correspond  more  closely  with  a seated  position.  Subsequently,  , 

: ' Si 

several  radiographs  of  seated  pilots  were  obtained  and  digitized  ^ 

so  that  the  configuration  could  be  determined.  The  base  point 

coordinates  in  the  y-z  plane  for  the  ad  hoc  configuration  and  one  of  | 

i 

the  radiograph  configurations  is  given  in  Table  5.  The  overall  j 

I 

length  of  the  thoracolumbar  columns  for  this  data  corresponds  closely  | j 

to  that,  of  an  aveiaye  iuale.  Clauser  (19C0)  reported  that  in  a cample 
o£  13  male  cadavers,  the  mean  distance  from  omphr-lian  to  cervicale 

I 

is  43.5  cm  and  cited  an  earlier  unpublished  Air  Force  study  in  which  j 

I 

the  same  figure  was  reported.  The  nearest  c'^r responding  dimension  in 
this  model,  the  distance  from  L4  to  Tl,  is  43.3  cm.,  which  is  in  good 
agreement.  { 

The  vertebral  body  heights  and  the  disc  heights  along  the  central 
axis  are  based  on  the  data  of  Lanier  and  Todd  and  Pyle  (1928) , in 
the  thoracolumbar  spine.  In  the  cervical  spine,  che  dimensions  are 
based  on  data  provided  by  AMRI,  and  the  orientations  were  measured  from 
radiographs. 
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Table  5 Geometry  of  Spine 


Vertebral 

Level 

Coordinates  of 
Inferior  End  Plate  Center 

y (cm)  z (cm) 

Vertebral  Body 
Height 

(cm) 

Intervertebral 
Disc  Height* 

(cm) 

L5 

1.800 

2.020 

2.392 

1.859  . 

L4 

1.100 

5.700 

2.636 

1.354 

L3 

1.000 

9.550 

2.751 

1.223 

L2 

1.331 

13,450 

2.792 

1.173 

LI 

2.142 

17.150 

2.726 

0.996 

T12 

3.003 

20.590 

2.567 

0.822 

Til 

3.882 

23.680 

2.433 

0.645 

TIO 

4.594 

26.500 

2.298 

0.477 

T9 

4.849 

29.240 

2.146 

C.460 

T8 

4.638 

31.830 

2.073 

0.459 

T7 

4.580 

34.300 

2.019 

0.404 

T6 

4.250 

36.610 

1.990 

0.314 

T5 

3.990 

38.850 

1.957 

0.266 

T4 

3.690 

41.000 

1.902 

0.214 

T3 

3.350 

43.150 

1.850 

0.274 

T2 

2.920 

45.260 

1.790 

0.306 

T1 

2.410 

47.440 

1.648 

0.448 

C7 

1.909 

49.420 

1.  612 

0.394 

C6 

1.760 

51.448 

1.516 

0.434 

C5 

1.460 

53.516 

1.515 

0.576 

C4 

1.290 

55.439 

1.513 

0.417 

C3 

1.484 

57.332 

1.511 

0.398 

C2 

1.636 

59.239 

1.500 

0.408 

* Indicates  disc  below  vertebral  level. 
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The  secondary  nodes  of  the  vertebrae  are  interconnected  by 
deformable  elements  which  represent  the  connective  tissues,  ligaments 
and  the  intervertebral  disc.  The  center  line  of  each  intervertebral 
disc  element  connects  the  centers  of  the  endplates  of  adjacent 
vertebra.  Each  intervertebral  disc  has  stiffness  in  bending  due  to 
both  flexion-extension  and  lateral  bending,  in  torsion,  and  in  the 
axial  mode. 

The  stiffnesses  are  based  primarily  on  the  experimental  work  of 
Markolf  (1970) , Brown,  Hansen  and  Yorra  (1957) , Rolander  (1966)  and 
Kazarian  (1972) . These  experimental  measurements  were  augmented  by 
model  studies  of  the  intervertebral  disc  under  axial  load  performed 
by  Kulak  (1974).  A large  part  of  the  procedure  of  assigning  these 
stiffnesses  has  been  summarized  in  Schultz,  Belytschko,  et  al  (1973). 
That  paper  also  lists  12  other  sources  of  intervertebral  disc  data 
which  have  been  considered  in  developing  the  daca  used  in  the  model. 

In  spite  of  the  wealth  of  literature  available  on  this  topic, 
the  experimental  data  is  inadequate  for  assigning  stiffnesses  at  each 
vertebral  level.  However,  geometrical  dimensions  at  every  level  have 
been  provided  by  Lanier  (1939) . In  order  to  estimate  stiffnesses 
at  each  level,  it  was  assumed  that  stiffnesses  vary  in  relation  to 
the  geometry.  This  is  equivalent  to  assuming  that  intrinsic  material 
properties  do  not  vary  significantly  with  vertebral  level.  Thus,  by 
using  the  experimental  values  at  the  levels  at  which  they  are  available 
and  by  using  the  geometrical  ratios  for  other  levels,  stiffnesses 
were  assigned  tc  every  level  of  the  spine. 

The  procedure  by  which  this  was  done  follows.  First,  each  disc 
was  idealized  as  an  elliptical  ring  of  linear,  elastic  material 
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corresponding  to  the  annulus  fibrosis.  The  outside  ellipse  is 
defined  by  major  and  minor  diameters  a^  and  b^,  respectively,  and 
the  inner  ellipse  by  major  and  minor  axes  a^  and  b^,  respectively. 
These  two  ellipses  were  assumed  to  be  concentric.  Lanier's  mean 
values  of  the  transverse  and  sagittal  diameters  of  each  vertebral 
body  were  assigned  to  a^  and  b^,  respectively,  and  to  the  height  h. 

The  interior  diameters  were  assumed  to  be  0.75  of  the  outer  diameters, 
based  on  Farfan  et  al.  (1970) . It  was  assumed  that  this  elliptical 
ring,  representing  the  annulus  fibrosis  and  longitudinal  ligaments, 
provides  the  major  resistance  to  bending,  torsion  and  shear.  This 
assumption  is  realistic  in  view  of  the  hydrostatic  behavior  of  the 
nucleus  pulposus  reported  by  Nachemson  (1960) , for  a hydrostatic 
material  does  not  contribute  to  shear  or  torsional  resistance,  and 
contributes  to  bending  resistance  only  if  it  is  located  asymmetrically 
with  respect  to  the  neutral  bending  plane.  The  axial  stiffness  was 
assumed  to  be  proportional  to  the  total  cross-sectional  area  of  the 
disc. 

It  follows  then  from  strength  of  materials  that  the  stiffnesses 
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vary  from  level  to  level  in  proportion  to  the 
factors:  for  tension,  compression,  and  shear 

following  geometric 

r,  = ^ (a  b - a.b.  ) 

A h o o 11 

for  bending  about  the  minor  diameter 

(4.1) 

r = ^ (a^b  - a?b . ) 
m n o o 11 

(4.2) 

87 

■THt  1 ■*’■' “M*  W|4», ^ ..,. 


for  bending  about  the  major  diameter 


1 , V.3  .3, 

r = r (a^b  - a.b. ) 

m h o o 11 


(4.3) 


for  torsion 


/ 3,  3 

= k( 

o o 


at+bT 
1 1 


(4.4) 


Extensive  data  is  available  on  the  stiffnesses  at  the  L3/L4  and 
T8/T9  levels  from  cadaver  material  measurements.  There  is  also  some 
data  on  the  upper  thoracic  spine.  From  the  data  the  disc  stiffnesses 
at  other  levels  were  calculated  by  using  the  geometrical  ratios 
listed  above.  Finally,  measured  stiffnesses  at  several  other  levels 
were  used  to  check  on  the  consistency  of  the  data,  which  was  found 
to  be  quite  good. 

Because  the  discs  are  very  short,  shear  action  is  quite  strong. 
To  account  for  this,  shear  factors  (|>  which  appear  in  Eq.  (1.4)  were 
added  to  account  for  the  shear  behavior. 

Another  aspect  which  had  to  be  considered  in  this  study  is  the 
nonlinearity  of  the  force  deflection  characteristic  of  the  disc.  Very 
little  data  is  available  from  which  the  nonlinearity  could  be  pre- 
cisely established,  except  for  perhaps  the  axial  force-deflection 
behavior  of  the  disc.  The  nonlinearity  in  the  axial  mode  is  distinctly 
different  in  the  lumbar  and  throacic  regions.  As  can  be  seen  from 
Fig.  12,  the  axial  behavior  of  the  disc  is  quite  nonlinear,  parti- 
cularly if  one  considers  both  the  tensile  and  compressive  regimes. 
However,  if  the  preloaded  state  resulting  from  body  weight  is 
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L+4000 

Figure  12.  Response  of  L3  to  axial  load. 


considered  as  a point  of  reference,  then  the  behavior  of  the  disc 
for  compressive  loads  above  that  point  is  rather  linear.  In  addition, 
it  can  be  seen  that  a linear  approximation  to  the  behavior  above  that 
point  requires  a greater  stiffness  coefficient  than  from  the  unloaded 
state,  by  a factor  of  about  1.33.  However,  the  data  used  in  the 
model  does  not  include  this  factor.  Although  the  model  could 
handle  a cubic  force  deflection  curve  which  would  approximate 
the  behavior  shown  in  Fig.  12  quite  closely,  it  was  not  used  because 
of  the  difficulties  of  introducing  the  preloaded  state  into  the  cubic 
curve . 

The  axial  force  deflection  characteristics  in  a thoracic  spine 
are  illustrated  in  Fig.  13.  It  can  be  seen  that  both  experimental 
measurements  and  computer  model  results  based  on  the  nonlinear  stress- 
strain  curve  as  reported  by  Kulak,  et  al.  exhibit  very  linear  behavior. 
Therefore  in  the  thoracic  and  cervical  regions,  a linear  curve  is 
adequate  and  no  further  modifications  are  contemplated.  It  may  be 
mentioned  that  the  increased  linearity  of  the  axial  behavior  of  the 
thoracic  and  cervical  discs  is  related  to  the  reduced  relative  height 
of  these  discs,  where  the  relative  height  is  defined  as  the  ratio  of 
height  to  average  diameter.  This  relative  height  is  considerably 
smaller  in  the  thoracic  and  cervical  regions.  For  discs  with  small 
relative  height,  the  application  of  axial  load  does  not  result  pri- 
marily in  hoop  stress  in  the  annulus,  but  a combination  of  compressive 
and  tensile  stresses,  so  that  the  nonlinearities  of  the  material  tends 
to  cancel.  The  damping  coefficients  were  based  on  the  work  of 
Payne  (1971)  and  Kazarian  (1972).  A summary  of  all  material  properties 
for  the  discs  is  given  in  Table  6. 
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: elements  were  axially  damped  using  a 0.2%  stiffness  proportiona 
damping. 


The  ligaments  are  represented  by  spring  elements  which  connect 
the  secondary  nodes  representing  the  spinous  process  tips  (inter 
and  intra  spinous  ligaments)  , spring  elements  connecting  the  transverse 
process  tips  (inter-transverse  ligaments)  and  springs  connecting 
points  on  the  vertebral  body,  which  represent  mainly  the  ligamenta 
flava.  All  of  these  elements  are  nonli'iear  in  that  they  have  no  resistance 
in  compression,  i.e.,  when  they  are  slack.  In  the  thoracolumbar  spine, 
the  articular  facets  are  also  represented  by  spring  elements.  The 
locations  of  the  secondary  nodes  to  which  these  spring  elements  are 
connected  were  adjusted  so  that  the  line  of  action  of  the  facet  forces 
corresponds  to  the  normal  to  the  facet  planes.  Typical  lines  of 
action  for  articular  facets  in  the  lumbar  and  thoracic  spine  are  given 
in  Table  7. 

The  model  of  the  cervical  spine  is  similar  to  that  of  the  thoraco- 
lumbar spine,  except  that  no  ligaments  were  included  at  this  time 
and  that  each  articular  facet  plane  was  represented  by  three  points, 
so  that  the  hydrodynamic  element  could  be  used  to  model  the  facet 
joint.  The  intervertebral  disc  heights  in  the  cervical  spine  were 
obtained  from  data  provided  by  AMRL.  The  data  included  both  an  anterior 
and  posterior  disc  height,  so  these  were  averaged  for  use  in  the  model. 
The  vertebral  geometry  was  obtained  from  direct  measurement  of  a 
cadaver  spine.  These  measurements  included  vertebral  body  height, 
endplate  areas,  spinous  process  location  and  three  points  on  each 
superior  articular  facet  plane.  Both  right  and  left  superior  articular 
facet  planes  were  measured  and  these  were  then  averaged  to  obtain  a 
symmetric  configuration.  The  inferior  articular  facet  planes  were 
obtained  by  placing  the  cervical  vertebrae  in  a standard  configuration 


Table  7.  Lines  of  Action  of  Articular  Facets 
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Vertebral  Level 

X- component 

y- component 

z-component 

L5-L4 

-0.96 

0.2 

0.18 

L4-L3 

-0.87 

0.48 

0. 

L3-L2 

-0.85 

0.42 

0.31 

L2-L1 

-0.99 

0.09 

0.08 

L1-T12 

-0.71 

0.54 

0.43 

T12-T11 

-1.00 

0. 

0. 

Tll-TlO 

0.13 

0.97 

0.14 

T10-T9 

0.13 

0.89 

0.44 

T9-T8 

0.16 

0.92 

0.33 

T8-T7 

0.14 

0.93 

0.31 

T7-T6 

0.39 

0.76 

0.51 

T6-T5 

0. 

0.99 

0.13 

T5-T4 

0. 

0.99 

0.07 

T4-T3 

0.44 

0.69 

0.56 

T3-T2 

0.22 

0.81 

0.53 

T2-T1 

0.31 

0.69 

0.65 
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obtained  from  a radiograph  and  then  defining  the  points  of  the 
inferior  articu.lar  facet  planes  so  that  they  are  parallel  and 
a prescribed  distance  from  the  superior  facet  planes;  this  dis- 
tance was  measured  in  a direction  normal  to  the  facet  plane.  The 
resulting  configuration  is  illustrated  in  Fig.  14.  The  facet 
planes  in  these  figures  appear  as  triangles.  The  force  deflection 
characteristics  of  the  facet  ioint  was  then  modelled  with  the 
hydrodynamic  element  described  in  Chapter  II. 


2.  Rib  Cage 

The  model  of  the  rib  cage  consists  of  rib  pairs  1-10  and  the 
sternum.  Each  rib  is  modelled  as  a rigid  body , with  the  deformability 
of  the  rib  cage  represented  by  deformable  elements  which  represent 
the  costo-transverse  and  costo-vertebral  articulations.  The  location 
of  the  superior  and  inferior  costo-vertebral  CV  articulations  and 
the  costo-transverse  CT  articulations  were  based  on  those  reported  in 
Schultz,  Benson  and  Hirsch  (1974).  The  diagram  from  Andriacchi,  etal. 
(1974)  indicates  the  nature  of  these  data.  The  geometry  of  each  rib 
is  characterized  by  8 points  placed  as  follows:  two  coincident 

points  in  the  costo-tuberco  defining  the  position  of  the  costo-trans- 
verse articulation,  a pair  of  inferior  and  superior  points  placed  at 
the  rib  head  at  the  pooitions  of  the  costo-vertebral  articulations, 
two  points  placed  along  the  inferior  and  superior  borders  of  the  rib 
shaft  at  the  mid  axillary  line  to  provide  points  of  attachment  for 
deformable  elements  representing  the  elastic  behavior  of  soft  tissues 
occupying  the  intercostal  spaces, and  a pair  of  points  at  the  interior 
end  of  the  calcified  portion  of  the  rib  providing  points  of  attachment 
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for  deformable  elements  which  represent  the  costo-cartilages.  Rib 
geometry  data  were  obtained  from  measurements  reported  in  Schultz, 
Benson  and  Hirsch,  (1974a).  Data  describing  the  rib  shaft  geometry 
are  reported  in  Table  9,  taken  from  Andriacchi,  et  al.  The  rib  cage 
geometry  is  illustrated  in  Figs.  8 and  9.  The  stiffness  properties 
of  the  deformable  elements  used  in  the  rib  cage  are  summarized  in 
Table  10. 


3.  Viscera  and  Abdominal  Cavity 

The  abdominal  cavity  and  viscera  are  represented  by  hydrodynamic 
elements  stacked  in  series  with  rigid  bodies  between  each  layer  as 
shown  in  Fig.  10.  The  compressibility  of  the  model  from  an  analytical 
viewpoint  is  thus  entirely  axial.  In  the  actual  viscera,  the  wave 
motion  is  mainly  governed  by  the  interaction  of  the  walls  of  the 
torso  and  its  contents.  Thus,  in  response  to  a compressive  load,  the 
viscera  would  move  vertically  and  laterally,  with  the  latter  compon- 
ent stretching  the  abdominal  walls. 

The  mechanism  of  wave  propagation  through  the  viscera  by  an 
interaction  of  the  membrane  lining  and  the  hydrodynamic  contents  has 
been  studied  by  Torvik  (1970) . He  derived  relationships  for  both 
large  deformations  of  the  membrane  and  nonlinear  membrane  response. 
For  small  deflection,  linear  membrane  response,  he  gives  the 
standard  water  hammer  formula  for  wave  speed 


.Jlil 

T 2ro 


where  E is  Young's  modulus  for  the  membrane  wall,  t the  thickness  of 
the  wall,  r its  radius  and  p the  density  of  the  fluid,  which  in  this 
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case  are  the  viscera.  This  formula  neglects  the  compressibility  of 

the  viscera.  If  we  let  E,  the  modulus  of  skin  (although  the 

6 2 

passive  resistance  of  muscles  is  also  included)  be 60  x 10  dynes/cm  , 
and  let  t=l  cm,  r=10  cm,  we  obtain  c = 1700  cm/sec.  This  corresponds 
very  well  with  the  values  measured  by  Weis  and  Mohr  (1967) , which 
are  1980  cm/sec. 

The  bulk  modulus  is  then  computed  by  the  standard  formula 
for  acoustic  wave  speed. 

B = pc^  (4.6) 

6 2 

This  gives  an  effective  value  for  B of  3.92  x 10  dynes/cm  . 

This  bulk  modulus  should  not  be  interpreted  as  a bulk  modulus  of 
the  viscera?  it  represents  a constant  that  reflects  the  combined 
action  of  the  membrane  walls  and  the  abdomen. 

The  visceral  elements  are  connected  through  rigid  bodies  to  ribs 
TIO.  This  implies  that  all  axial  load  in  the  abdominal  cavity  is 
transferred  to  the  rib  cage;  no  axial  load  transfer  to  the  interior 
of  the  thoracic  cavity  is  assumed.  The  bottom  visceral  elements  are 
connected  to  the  pelvis. 


4 . Head 

The  head  is  modelled  as  a rigid  body.  In  the  examples  studied, 
any  helmet  or  helmet  mounted  devices  were  assumed  rigidly  attached 
to  the  head,  so  the  mass  and  mass  moments  of  inertia  of  these  were 
simply  added  to  the  head.  However,  it  is  possible  to  separate  the 
inertias  of  the  head  and  helmet  and  interconnect  the  two  by  springs 
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as  shown  in  Fig.  15. 

The  mass  of  the  head  was  chosen  to  be  4.6  kg,  the  mass  moments 

2 

of  inertia  about  400  kg-cm  based  on  the  anthropometric  dummy  mea- 
surements of  Bartz  (1972) . The  helmet  masses  were  based  on  the 

average  of  3 helmets  measured  by  the  investigators.  These  had  an 

2 

average  mass  of  1.4  kg  and  mass  moments  of  100  kg-cm  . 


5.  Pelvis 

The  pelvis  was  represented  as  a rigid  body  with  a mass  and  mass 
moments  of  inertia  given  in  Table  11.  The  data  is  taken  from  Bartz 
(1972) . The  geometrical  aspects  of  the  pelvic  representation  are 
shown  in  Figs.  8 and  9. 


6.  Preliminary  Evaluation  of  Iniury  Potential 


I .’^1 

I ' 


In  order  to  obtain  a qualitative  estimate  of  the  injury  potential 

of  various  combinr.tions  of  sagittal  plane  moments  and  axial  force 

for  the  vertebrae,  a method  of  calculating  the  maximum  stress  in  the 

vertebrae  based  o:.  the  combined  axial  force  and  moment  predicted  by 

the  model  was  developed.  This  is  only  a simplified  model,  but  it  is 

indicative  of  the  effects  of  moments  on  stress  levels  in  the  vertebral 

bodies.  Since  injury  is  probably  related  to  stress  levels,  this 

method  gives  a means  for  evaluating  the  effects  of  the  moments.  We 

idealize  the  vertebral  body  as  a cylindrical  shell  of  radius  r^  and 

height  h.  The  shell  is  considered  to  be  cortical  bone  with  a Young's 

11  2 

modulus  of  1.5  X 10  dynes/cm  . The  interior  of  the  shell  (vertebral 
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Head 


core)  is  trabecular  bone  (soft  bone) , with  a mean  modulus  of  elast- 

8 2 

icity  = 7.35  X 10  dynes/cm  These  moduli  were  obtained  from 
Evans  (1970).  The  following  aspects  of  the  vertebra's  geometry  have 
been  neglected: 

1.  The  variation  of  cross  section  with  z,  the  cross  section  of  the 
vertebral  body  as  shown  in  Fig.  16  is  simplified  as  shown. 

2.  The  ellipticity  of  the  endplates  is  neglected  and  a circular 
cross  section  is  assumed. 

The  maximum  stress  for  any  combination  of  moment  M and  axial 
force  P occurs  in  the  cortical  bone  and  is  obtained  by  superimposing 
the  stress  due  to  bending  Og  and  the  axial  stress  a^.  The  bending 
stress  is  computed  from 

o_=jE(r)  (4.7) 

tj  p 

where  p is  the  radius  of  curvature,  E is  Young's  modulus,  which  is 
given  by 


E(r) 


jEi  r < r. 
^ " ^i 


(4.8) 


ro  /-7r/2 


M 


= / ' r°L 


rcos0 


E(r)  rcos0  rdrd0 


o o 


^ fE.r.^  + E (r^  - r^)  } 
4p  ■ 1 1 o o i'  ^ 


(4.9) 


4Mr_E 

_ o o 

°max  TT{E.r.4  + r.'*)} 

11  O O 1 


+ 


TT{E.r . ^ 

1 1 


PEo 

+ E (r  ^-r . 2) } 
o o 1 


(4.16) 


An  equivalent  circular  radius  r^  for  the  vertebral  bodies  was 
determined  so  that  the  area  of  the  two  geometries  was  the  same.  The 
thickness  of  the  cortical  bone  was  assumed  to  be  a constant  0.3  mm 
for  all  the  vertebrae,  i.e.  maximum  axial  force 

for  each  vertebrae  were  obtained  from  Payne  (1971) , who  summarized 
the  results  of  several  experiments  ccncerned  with  the  compressive 
breaking  load  of  individual  vertebrae.  From  this  data,  the  maximum 
breaking  stress  under  only  axial  load  and  the  pure  moment  without 
axial  force  were  determined.  The  results  of  these  calculations  are 
summarized  in  Table  12.  As  can  be  seen,  the  breaking  strength  as 
ccxnputed  by  this  formula  is  fairly  constant  over  the  entire  spine, 
except  for  moderate  deviations  in  the  lumbar  and  upper  thoracic  verte- 
Figure  25  depicts  the  rooment”axial  force  interaction  line  of 
several  vertebral  levels  (alternating  levels  wei  i omitted  for  clarity). 
Each  point  on  a line  represents  the  combination  of  moment  and  axial 
force  which  corresponds  to  the  maximum  stress  for  that  vertebral 
level. 
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Inertial  properties  for  the  motion  segments  of  the  isolated 
thoracolumbar  spine  were  taken  from  Liu,  et  al  (1973)  and  are  given 
in  Table  11.  Inertial  properties  for  the  rib  cage  model  were 
assigned  by  attributing  10%  of  the  thoracic  mass  to  this  portion 
of  skeletal  structure.  For  the  viscera  model  the  lumbar  mass  was 
equally  distributed  among  the  five  vertebrae  and  two  rigid  bodies 
used  in  the  model.  The  sagittal  plane  rotational  inertia  of  the 
lumbar  vertebrae  were  reduced  based  on  measurements  and  calcula- 
tions obtained  from  anatomical  cross  sectional  geomer.ries.  A 
summary  of  the  inertial  properties  of  the  spine  torso  model  with 
rib  cage  is  given  in  Table  13.  Inertial  properties  for  the  cer- 
vical motion  segments  were  obtained  by  distributing  the  mass  uni- 
formily  at  each  vertebral  level;  see  Table  13. 
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Table  15.  Inertial  Properties  for  Complete  Model 


Motion 

Segment 

Mass  Pj 
Grams  x 10^ 

1 , 
xxl 

Gram-cm^  x 10^ 

1 , 

5 

Gram-cm'^  x 10 

^zzl 

Gram-cm^  x ; 

Pelvis 

16.200 

128.000 

20.000 

19.300 

L5 

1.500 

2.783 

1.795 

2.382 

L4 

1 . 500 

2.748 

1.704 

2.291 

L3 

1.500 

2.809 

1.682 

2.280 

L2 

1.500 

2.840 

1.695 

2.291 

LI 

1.500 

2.740 

1.569 

2.212 

T12 

1.556 

7.002 

1.309 

1.919 

Til 

1.453 

7.056 

1.230 

1.941 

TIO 

1.202 

6.028 

1.129 

1.648 

T9 

1.267 

6.164 

1.230 

1.716 

T8 

1.176 

5 543 

1.208 

1.670 

T7 

1.158 

S.347 

1.219 

1.659 

T6 

1.043 

4.425 

1.162 

1.546 

'i  5 

1.025 

3.383 

1.151 

1.490 

T4 

0.964 

3. 138 

1.060 

1.354 

T3 

1.010 

2.878 

1.174 

1.422 

T2 

0.974 

2.007 

1.029 

1.230 

T1 

1.209 

0.745 

0.518 

1.716 

C7-C2 

1.000 

0.700 

0.500 

1.500 

Head 

5.612 

44.786 

4.044 

3.585 

Ribs  T1 
to  TIO 

0.074 

0.573 

0.074 

0.074 

Lower 

Viscera 

1.500 

10.700 

0.550 

1.000 

Upper 

Viscera 

1.500 

10.700 

0.550 

1.000 

11 


8.  Summary  of  Models 

Model  I,  shown  in  Figure  17,  consists  of 

i.  thoracolumbar  spine 

ii.  a single  beam  element  which  represents  the 
cervical  spine 

iii.  pelvis 

iv.  head 

Model  II,  shown  in  Figure  18,  consists  of 

i.  thoracolumbar  spine 

ii.  a single  beam  element  for  cervical  spine 

iii.  rib  cage 

iv.  viscera  represented  by  hydrodynamic  elements 
V.  pelvis 

vi.  head 

Model  III,  shown  in  Figure  19,  consists  of  ' 

j 

i.  thoracolumbar  spine 

ii.  cervical  spine  modelled  as  individual  vertebrae  | 

iii.  pelvis 

iv.  head 
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CHAPTER  V 


STUDIES  OB’  SPINAL  RESPONSE 


1.  Isolated  Ligamentous  Spine  Model 


The  first  group  of  studies  reported  here  were  conducted  with 
the  ^odel  of  the  isolated  ligamentous  spine  described  as  Model  I 
in  the  previous  chapter.  All  of  these  problems  were  studied  as 
three  dimensional  problems,  even  though  in  most  cases  a two  dim- 
ensional model  would  have  been  adequate.  In  all  cases,  explicit 

-4 

integration  was  used  with  a time  increment  of  10  seconds,  re- 
quiring 800  time  steps  for  a solution  at  80  msec. 


Response  of  unrestrained  spine  under  10  The  first  study  was 

made  to  evaluate  the  response  of  the  isolated  ligamentous  spine 
under  a perfectly  vertical  lOG^  acceleration.  A rate  of  onset  of 
714  g/sec  was  used;  the  acceleration  was  then  maintained  at  a con- 
stant level  of  lOg  for  66  msec.  Unlike  other  studies  described  in 
this  section,  no  seatback  or  restraints  were  included  in  this 
solution. 

The  deformed  configurations  at  40,  60  and  80  msec  are  shown 
in  Figs.  20,  21,  and  22.  As  can  be  seen,  particularly  at  80  msec, 
severe  curvatures  of  the  spine  have  developed,  and  the  spinal 
column  can  be  considered  to  have  buckled.  Once  the  column  has 
buckled,  the  computed  response  is  quite  unrealistic.  For  example, 
moments  on  the  order  of  3000  N-cm  are  developed  in  the  lumbar 
spine. 

Liu,  et  al  (1973)  have  observed  a similar  response  in  a homo- 
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geneous  bar  model  of  the  spine,  but  were  inconclusive  as  to 

whether  this  response  is  correct.  However,  if  one  considers  that 

Lucas  and  Bressler  (1961)  have  shown  a static  frontal  plane 

buckling  load  of  20  to  100  Newtons  for  the  thoracolumbar  spine 

constrained  against  sagittal  plane  buckling,  then  this  response 

is  not  unexpected.  The  bending  stiffness  of  the  isolated  spine  in 

the  sagittal  plane  is  of  the  same  order  as  the  bending  stiffness 

in  the  frontal  plane,  and  the  initial  sagittal  plane  curvature 

of  the  spine  reduces  its  sagittal  mode  buckling  even  further.  The 

lOG  environment  results  in  axial  loads  of  4500  Newtons,  and  the 
z 

duration  is  on  the  order  of  200  msec,  which,  as  shown  subsequently, 
is  of  the  order  of  the  lowest  flexural  periods  and  4 times  the 
axial  period  of  the  spinal  column,  so  inertial  effects  cannot  pre- 
clude buckling.  Hence  it  is  not  unexpected  that  the  unrestrained 
spine  will  buckle  in  this  environment. 

There  are  essentially  two  mechanisms  that  prevent  buckling 
in  the  actual  spine: 

i.  the  action  of  the  restraint  system,  seacj'ack  and 
musculature; 

ii.  the  interaction  with  the  torso  and  rib  cage,  which  may 
significantly  increase  the  bending  stiffness  of  the  composite 
spine/torso. 

The  results  ir.  subsequent  Sections  reflect  rome  attempts  to 
study  the  effects  of  the  second  mechanism.  In  all  studies  describe' 
in  the  remainder  of  this  section,  the  seatback,  restraint  system, 
and  a simplified  torso  representation  were  included.  These  were 
sufficient  to  eliminate  buckling. 
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An  alternative  to  including  constraints  in  the  model  would  be 
to  use  a small  displacement  method  of  analysis.  However,  small 
displacement  analyses  are  valid  only  as  long  as  the  actual  displace 
ments  of  the  spine  are  small.  In  problems  such  as  that  of 
eccentric  head  loading,  large  displacements  are  undoubtedly  encount 
ered  and  would  invalidate  a small  displacement  analysis. 

Another  point  of  interest  in  the  response  of  this  model  is 
that  the  effects  of  the  ligaments  are  minimal.  Because  of  the 
large  compression  of  the  discs,  the  ligaments  as  modelled  here  all 
become  slack  and  provide  no  stiffness.  This  situation  continues 
even  as  buckling  is  initiated:  the  ligaments  do  ,iOt  prevent 
buckling.  This  is  also  true  of  the  spring  models  of  the  articular 
facet  joint. 

Effect  of  sagittal  plane  curvature.  Two  distinct  curves 
characterize  the  sagittal  plane  curvature;  the  lumbar  curve  formed 
by  the  vertebrae  L5  through  LI,  and  the  thoracic  curve  consisting 
of  T12  through  Tl.  Sagittal  plane  curvature  changes  from  concave 
!,  j convex  (when  viewed  posteriorly)  in  the  region  between  L2  and 
Til.  In  this  study  three  magnitudes  of  sagittal  curvature  were 
simulated.  The  first,  shown  in  Fig.  23,  represents  a seated  con- 
figuration in  which  the  lumbar  curvature  averages  0.045  cm  ^ and 
the  thoracic  curvature  0.017  cm  This  configuration  was  chosen 
to  illustrate  an  occupant  with  large  lumbar  curvature.  The 
second,  shown  in  Fig.  24,  is  an  erect  configuration  in  which  the 
average  lumbar  curvature  is  0.025  cm  This  configuration  is 
based  on  several  radiographs  of  pilots  in  ejection  seats.  The 
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third,  shown  in  Fig.  25,  represents  a hypothetical  configuration 
in  which  all  sagittal  plane  curvature  has  been  eliminated.  In 
this  "straight  spine"  model  all  of  the  intervertebral  disc  ele- 
ments are  along  the  vertical  axis  as  are  the  mass  centers. 

The  maximum  axial  forces.  Fig.  26,  for  the  three  configura- 
tions decrease  smoothly  from  bottom  to  top,  particularly  in  the 
"straight  spine"  model,  and  the  curvature  has  little  effect  on  the 
values.  On  the  other  hand,  a comparison  of  maximum  sagittal  plane 
moments  for  the  three  models.  Fig.  27,  shows  large  differences 
among  the  models.  The  maximum  sagittal  plane  moment  in  both 
curved  spines  occurs  at  L1-T12,  the  point  where  spinal  curvature 
changes  from  concave  to  convex.  In  general,  the  sagittal  plane 
moments  for  the  erect  configuration  are  smaller  than  for  the 
curved  configuration  at  all  levels,  while  in  the  straight  spine, 
the  moments  are  smallest.  These  results  indicate  that  the  pre- 
ejection configuration  of  the  pilot  is  an  important  factor. 

While  the  compressive  forces  in  the  vertebrae  are  not  effected  by 
configuration  changes,  the  additional  stress  due  to  the  large 
bending  moments  acting  on  the  vertebrae  may  increase  the  possibil- 
ity of  injury. 

Figures  28  and  29  show  the  maximum  moment  and  axial  force 
predicted  by  the  model  for  both  of  the  curved  configurations  in 
the  vertebral  bodies,  to  illustrate  the  simplified  injury  poten- 
tial model  developed  in  Chapter  IV.  The  axial  force  and  moment 
values  shown  represent  the  average  of  the  axial  forces  and  moments 
for  each  vertebrae  in  the  discs  immediately  above  and  below  it. 
These  values,  in  all  cases,  are  the  maximum  that  were  observed 
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INTERVERTEBP 


MAXIMUM  AXIAL  FORCE 


A CURVED 
O ERECT 


Figure  26.  Comparison  of  maximum  axial  forces  in  spines  with 
and  without  sagittal  plane  curvature. 
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during  the  80  msec  simulation  (the  maocimum  axial  forces  and 
moments  do  not  necessarily  occur  at  the  same  time) . The  length  of 
the  line  segment  with  the  arrow,  relative  to  the  distance  of  the 
limit  load  from  the  origin,  is  one  indication  of  the  probability 
of:  no  injury  at  a level  (if  the  point  falls  inside  the  limit  line)  . 
Thus,  the  probability  of  injury  is  somewhat  greater  for  the  more 
curved  spine,  particularly  in  the  T1-T5  region,  where  the  distance 
to  the  injury  line  is  almost  20%  less  for  the  more  curved  spine. 
Particularly  if  we  consider  the  possibility  that  the  moments  may 
serve  as  a triggering  factor  for  failures,  then  the  curved  spine, 
with  its  greater  moments  appears  to  have  a greater  injury  poten- 
tial. 


Rate  of  Onset.  One  of  the  parameters  of  interest  in  the  ejection 
problem  is  the  effect  of  changes  in  the  rate  of  onset  for  the 
acceleration.  Rapid  rates  of  onset  result  in  dynamic  magnifica- 
tion of  force  and  acceleration  magnitudes.  Hess  and  Lombard  (1957) 
reported  that  onsets  of  about  60  msec  produced  considerably  higher 
acceleration  levels  in  the  body  than  the  maximum  acceleration  of 
he  ejection  seat.  Alternatively,  slower  onset  rates,  (100  msec) 
do  not  magnify  acceleration  levels,  but  may  not  provide  adequate 
displacement  of  the  pilot  to  clear  the  aircraft  in  the  required 
time.  Two  r£imp  acceleration  profiles  were  simulated  here.  In  both 
the  maximum  acceleration  is  lOg;  the  slower  onset  profile  reached 
the  maximum  acceleration  at  40  msec  with  an  onset  rate  of  250g/sec, 
the  other  profile  reached  maxx..'Mm  acceleration  at  14  msec  with  an 
onset  rate  of  714g/sec. 


m 


Comparison  of  the  maximum  acceleration  of  the  head  showed  a 
19.07g  peak  at  60  msec  for  the  rapid  onset  of  18.03g  peak  at  72 
msec  for  the  slow  onset,  which  is  a 5%  reduction  in  the  accelera- 
tion level.  Similarily,  a comparison  of  maximum  forces  in  the 
lumbar  region  showed  the  axial  disc  forces  were  reduced  by  about 
8%.  Also  the  sagittal  plane  moments  and  facet  forces  were  re- 
duced by  10%  to  15%.  Thus  the  mcdel  predicts,  as  e-.^pected,  that 
slower  rates  on  onset  reduce  the  dynamic  magnification  of  axial 
loads  and  accelerations  and  reduce  the  overall  bending  response 
of  the  spine.  Time  histories  of  the  head  acceleration  and  forces 
are  shown  in  Figs.  30  to  33. 

Angled  Pulse  (slanted  seat) . The  orientation  of  the  ejection  seat 
with  respect  to  the  direction  of  the  acceleration  vector  is  an- 
other consideration  in  the  ejection  problem.  By  reclining  the 
seat  so  that  the  acceleration  vector  has  a slight  anterior  compon- 
ent with  respect  to  the  axis  of  the  spine,  two  beneficial  effects 
arc  introduced:  (1)  the  acceleration  component  along  the  axis  of 

the  spine  is  reduced,  and  (2)  the  resulting  anterior  component  of 
acceleration  provides  support  for  the  spine  by  forcing  the  seat- 
back  against  the  torso.  Two  orientations  of  the  acce Iteration 
vector  were  studied,  one  in  which  the  axis  of  the  spine  and  accel- 
eration vector  were  aligned  vertically  with  a peak  acceleration 
magnitude  of  lOg,  and  the  other  with  a 30°  included  angle  and  a 
lOg  peak  acceleration  magnitude  resulting  in  8.66g  component  along 
the  axis  of  the  spine  and  a 5g  anterior  component. 

Table  14  compares  the  maximum  lumbar  forces  for  three  simula- 
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Figure  33.  Comparison  of  moment  in  T4/T3  for  fast  and  slow  rate 
of  onset. 
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tions  of  80  msec  duration:  rapid  onset  with  30°  included  angle, 
slow  onset  with  no  angle,  and  rapid  onset  with  no  included  angle. 
These  results  show  a similarity  in  response  for  the  rapid  onset 
at  a 30°  angle  and  the  slow  onset  at  a 0°  angle,  while  the  rapid 
onset  at  0°  exhibits  generally  higher  force  levels.  The  similarity 
in  results  for  the  rapid  onset  at  30°  and  the  slow  onset  at  0° 
can  be  attributed  to  a reduction  in  the  effective  inertial  forces 
in  booh  simulations.  In  the  slow  onset  simulation,  the  effective 
inertial  force  is  reduced  by  the  decrease  in  dynamic  magnification. 
In  the  rapid  onset  at  30°,  a similar  reduction  in  inertial  force 
results  from  the  action  of  the  seatback,  which  supports  part  of 
the  inertial  load. 

These  results  appear  to  indicate  that  within  the  constraints 
of  maximum  possible  seat  angle  and  the  required  displacement  of  the 
ejection  seat  for  clearance  of  the  aircraft,  an  optimum  combination 
of  seat  angle  and  rate  of  onset  could  be  determined.  This  optimum 
solution  should  reduce  force  levels  in  the  spine. 


Eccentric  Head  Loading.  An  eccentric  head  mass,  which  represents 
a helmet  mounted  device,  was  chosen  to  illustrate  the  behavior  and 
applicability  of  the  model  to  this  class  of  problems.  In  this 
study,  the  head,  in  addition  to  its  own  mass,  included  a rigidly 
connected  0.9  kg  mass  located  0.102  meters  from  the  sagittal  plane. 
Because  of  the  positioning  of  the  added  mass,  the  problem  is  not 
symmetric  about  the  sagittal  plane,  and  a three  dimensional 
analysis  is  needed. 

Figure  34  shows  the  axial  force  time  histories  at  three  levels 
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TIME  (milliseconds) 

Figure  34.  Response  of  spine  to  eccentric  mass  distribution 
of  head. 
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of  the  spine.  As  can  be  seen  from  Fig.  34  the  axial  behavior 
clearly  demonstrates  t!ie  progression  of  a compressive  wave  up  the 
spine  at  about  SOm/sec.  This  value  of  axial  wave  speed  in  the 
spine  is  in  good  agreement  with  the  values  reported  by  Hess  and 
Lombsrd  (1957)  and  Li,  et  al  (1970).  The  peak  compressive  forces 
occur  earlier  in  the  upper  levels  of  the  spine,  due  to  the 
reflected  expansion  wave  which  moves  down  from  the  head  cancelling, 
in  pare,  the  upward  moving  compressive  wave.  Figure  34  shows  the 
cervical  moment  in  the  frontal  plane  due  to  the  eccentric  loading. 
This  moment  vanishes  until  the  compressive  wave  has  reflected  from 
the  head;  after  reflection,  a frontal  plane  moment  is  generated 
in  the  cervical  spine.  This  moment  attains  its  maximum  at  about 
65  msec,  which  is  considerably  later  than  when  the  peaks  in  com- 
pressive forces  are  reached. 

In  the  cervical  spine,  the  intervertebral  discs  have  very 
little  bending  stiffness,  whereas  the  articular  facets  have  large 
moment  arms  about  the  sagittal  plane.  Hence,  the  forces  in  the 
facets  can  be  estimated  neglecting  the  intervertebral  disc  moments 
by  ascribing  the  moment  to  a difference  of  the  vertical  forces  in 
a facet  pair.  This  calculation  indicates  that  the  peak,  compres- 
sive facet  force  is  on  the  saime  side  as  the  eccentric  mass  and 

3 

that  its  magnitude  is  about  10  Newtons,  which  is  on  the  order  of 
two  to  three  times  the  facet  force  without  an  eccentric  head  mass. 
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The  studies  reported  in  this  Section  were  conducted  with  the 
I more  complex  model  of  the  spine  described  as  Model  II  in  the  pre- 


vious chapter.  With  the  exception  of  the  last  study,  the  seat- 


back  and  restraints  were  included  in  all  studies  and  the  models 

were  driven  by  prescribing  the  acceleration  of  the  ejection  seat 

to  be  with  a rate  of  onset  of  714g/sec  over  the  first  14 

msec  to  the  maximum  acceleration.  The  large  bandwidth  of  this 

model  requires  that  the  explicit  integration  technique  be  used.  A 

-4 

time  increment  of  10  was  sufficient  for  stability  in  energy. 


Response  of  Rib  Cage  Model  Under  10G„.  In  the  first  study,  the 

response  of  the  rib  cage  model  without  any  viscera  was  studied 

under  a vertical  lOG  acceleration.  Deformed  configurations  at  i 

i 

20,  40,  60,  and  80  msec  are  shown  in  Fig.  35.  As  can  be  seen,  the  j 

i 

entire  rib  cage  collapses  around  the  spinal  column.  There  are 

essentially  two  mechanisms  that  prevent  the  collapse  of  the  rib 

cage:  j 

i.  the  action  of  the  musculature; 

ii.  the  interaction  of  the  viscera  and  rib  cage;  as  the 
viscera  are  compressed  thej  subject  vertical  forces  on  the  dia- 
phragm which  supports  the  rib  cage. 

Although  no  attempt  was  made  to  represent  the  complex  inter- 
actions of  the  musculature,  consideration  of  the  second  mechanism 
led  to  the  development  of  the  previously  described  multi-elemant 
viscera  model.  As  the  results  in  the  following  section  show,  by 
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Figure  35.  Response  of  complex  model  without  viscera 


including  a representation  of  the  viscera  the  collapse  of  the  rib 


cage  was  prevented. 


Response  of  Complete  Model  Under  10G_.  This  study  was  made  to 

z 

determine  the  response  of  the  combined  rib  cage  and  viscera  models 
under  +G^  acceleration.  Figure  36  depicts  the  deformed  configura- 
tions of  the  model  at  20,  40,  60,  and  80  msec.  This  figure  also 
illustrates  that  the  collapse  of  the  rib  cage  described  previously 
is  prevented  by  incorporating  the  viscera  into  the  model. 

More  importantly,  spinal  response  is  altered  by  the  inclusion 
of  the  rib  cage  and  viscera.  A comparison  of  maximum  lumbar  forces 
for  the  complete  model  and  the  isolated  spine  (Table  15) , shows 
an  overall  reduction  in  force  levels.  Experiments  by  Tennyson  and 
King  (1974) , in  which  eviscerated  cadavers  were  subjected  to  +G^ 
acceleration  both  with  and  without  the  abdominal  cavity  pressurized 
show  a reduction  in  axial  force  levels  in  the  disc  and  articular 
facets  for  the  lumbar  region  when  tip  abdominal  cavity  was  pres- 
surized. In  particular,  the  axial  force  levels  are  reduced  by 
about  12%  to  23%  which  compares  with  the  10%  to  25%  reduction  re- 
ported in  the  experiments.  The  sensitivity  of  the  axial  force 
reductions  to  the  value  of  the  bulk  modulus  of  the  viscera  elements 

was  also  considered;  decreasing  the  modulus  by  a factor  of  four 
6 2 

(1x10  dynes/cm  ) resulted  in  a 5%  to  12%  reduction  in  force  levels 

7 2 

while  increasing  the  modulus  by  a factor  of  four  (1.6x10  dynes/an  ) 
reduced  the  axial  forces  by  20%  to  35%. 

The  maximum  internal  pressures  of  the  viscera  model  ranged 

5 2 5 2 

from  5.8x10  dyne/cm  in  the  lower  elements  to  2.5x10  dyne/cm  for 


EJECTION  SIMULATION  VilTH  PVP  VISCERA  MODEL 


Figure  36,  Response  of  complete  model  to  lOG 
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the  upper  visceral  elements.  These  values  of  internal  pressure 

5 2 

are  substantially  higher  than  the  1.4x10  dyne/cm  value. used  for 

pressurization  in  the  Tennyson  and  King  experiments , or  the  pr-pc- 
5 2 

sures  of  2x10  dyne/cm  reported  by  Morris,  et  al  (1961)  in 

static  measurements  during  weight  lifting  experiments.  However, 

it  is  felt  that  this  range  of  visceral  pressure  is  not  unrealistic. 

Using  the  counterpart  of  Eg.  (4.5)  and  the  data  employed  to  obtain 

this  wave  speed,  the  maximxim  stresses  in  the  abdominal  wall  are 

6 2 

computed  to  be  5.8x10  dyne/cm  (84  psi) . These  stresses  could  be 

confined  by  the  musculature  and  the  tissues  alone,  for  the  strength 

8 2 

of  these  tissues  is  on  the  order  of  10  dyne/cm  . The  resulting 
expansion  of  the  cavity  would  be  0.1  cm,  which  is  probably  unde- 
tectable. 

The  addition  of  the  torso  to  the  model  reduced  the  bending 
moments  in  the  thoracic  spine.  Thus,  although  the  abdominal 
cavity-viscera  model  has  almost  no  inherent  bending  stiffness, 
the  overall  effect  of  this  additional  c ilumn  is  to  stiffen  the 
model  in  bending.  This  behavior  is  also  reflected  in  the  accelera- 
tion of  the  head.  The  maximum  acceleration  is  19. 6g,  which  is 
slightly  higher  than  that  of  the  isolated  column  for  the  same  con- 
ditions . 

Response  of  the  Rib  Cage  Model  During  Frontal  Impact.  To  illus- 
trate another  application  for  the  dynamic  response  spine  model,  a 
simulation  of  frontal  impact,  as  would  occur  in  an  automobile 
impacting  a barrier,  was  conducted.  The  spinal  model  including 
head,  pelvis,  legs  and  rib  cage  was  given  an  initial  velocity  of 
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V-:; 
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1.34x10  cm/sec  (30  inph)  in  the  anterior  direction,  with  the  pelvis 

4 2 

subjected  to  a 6.7x10  cm/soc  (68. 3g)  posterior  deceleration  for 
20  msec.  The  deceleration  pulse  was  a step  function  with  a de- 
creasing ramp  at  the  end  of  the  pulse,  the  magnitude  of  the  decel- 
eration corresponds  to  an  experimental  measurement  of  a 30  mph 
crash. 

Displacement  of  the  pelvis  was  prescribed  to  be  consistent 
T?ith  the  initial  velocity  and  deceleration,  such  that  the  displace- 


ment from  the  initiation  of  the  deceleration  to  20  msec  was  13.42 

I 

I 

cm.  The  pelvis  was  constrained  from  rotation  in  the  sagittal  | 

plane  to  simulate  the  effect  of  a lap  belt  restraint,  although  the  | 

elasticity  of  the  restraint  belt  was  not  included  in  this  Simula-  j 

tion . 

The  head  displaced  52.85  cm  anteriorly  at  40  msec  and  was 
4 2 

subjected  to  1.63x10  cm/sec  (16.65g)  peak  acceleration,  also  at 
40  msec.  An  axial  expansion  wave  travels  up  the  spine  at  about 
30  m/sec  which  is  the  same  wave  speed  calculated  for  the  com- 
pressive wave  in  the  ejection  simulations.  Peak  tensile  axial 

0 

forces  in  the  intervertebral  disc  ranged  from  b. 84x10  dynes  at 

8 8 
the  sacrum  L5  level  to  7.82x10  dynes  at  Tll-TlO  to  4.165x10 

dynes  at  the  T3-T2  level.  Figure  37  shows  the  spinal  configura- 

i 

tions  initially  and  at  20  and  40  msec.  The  exaggerated  position  j 

of  the  upper  torso  at  40  msec  is  due  to  the  fact  that  no  shoulder  j 

j 

harness  restraint  was  included  in  the  simulation  and  the  possible  i 

interaction  of  the  thoracic  cavity  with  the  steering  column  was  1 

! 

I 

not  considered.  It  also  appears  that  the  stiffness  of  the  ribs  1 

against  rotation  is  too  low.  | 
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Figure  37.  Frontal  impact  deformed  configurations 


PH  |I  ijHf  I 


3.  Modal  Analyses 

In  ordei  to  gain  a better  understanding  of  the  dynamic 
response  properties  of  the  spine,  modal  analyses  were  performed 
on  the  ligcunentous  isolated  column.  The  first  set  of  results  in 
Table  16  gives  the  lowest  7 natural  frequencies  of  the  isolated 
ligamentous  spine.  Model  I.  In  the  modal  analyses,  neither  the 
seat  nor  the  harness  restraints  are  included.  The  spine  was  con- 
strained from  motion  in  the  frontal  plane;  only  motions  in  the 
sagittal  plane  were  considered.  Material  properties  that  do  not 
account  for  the  preload  of  the  body  were  used. 

As  can  be  seen  from  the  results,  the  model  is  characterized 
by  a large  number  of  very  low  natural  frequencies.  All  of  the 
natural  modes  associated  with  these  frequencies  involve  primarily 
bending  deformations  of  the  spine.  The  lowest  natural  mode  of 
the  system  with  significant  axial  deformation  is  the  7th  mode, 
which  has  a frequency  of  17  cps. 

When  the  harness  was  included,  the  natural  frequencies  shifted 
as  shown.  Adding  the  harness  introduces  additional  frequencies 
which  correspond  to  the  motion  of  the  pilot  relative  to  the  seat. 
The  lowest  such  mode  is  vertical,  for  the  harness  exerts  little 
vertical  constraint.  The  fourth  frequency,  5.62  cps,  agrees  rea- 
sonably with  resonance  peaks  found  in  the  driving  seat  impedance 
measurements  by  Vogt,  et  al  (1968) . However,  the  bending  frequen- 
cies seem  to  be  sensitive  to  the  nature  of  the  model,  so  it  is  not 
clear  that  this  correlation  is  definitive. 

In  order  to  investigate  these  results  further,  the  following 
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Table  16- 


M^<-nral  Frequerun^s__oL-S£Hl£j^^^^ 


Without 

Harness 


1 

2 

3 

4 

5 

6 
7 


1.28 

3.14 

5.99 

9.94 

13.31 

16.71 

18.45 


Model  I 

With 

Harness 


c.q^htal  Plane  Freauencies 

Prasad-King  Model 


0.31 

1.78 

2.35 

5.62 

10.25 

13.08 

14.73 


Without 
Head  or  Hips 


3.66 

8.17 

13.01 

17.70 

22.21 

26.92 

30.90 


1.38 

4.48 

7.70 

11.74 

14.69 

23.30 

26.45 


1 

2 

3 

4 

5 

6 


av-ial  Frequenci^ 

Prasad-King  Model 


W i th  Ileaci 
and  Hips 


Model  I 

Without 
Head  or  Hips 


17.09 

32.31 

51.29 

77.22 

100.74 

124.89 


30.12 

19.41 

61.02 

39.29 

87.08 

80.61 

114.87 

113.18 

142.42 

148.81 

167.41 

179.48 

150 


yrr^.^  . 
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additional  modal  anal^^ses  were  made 

i.  axial  analyses  in  which  the  model  was  constrained  from 
all  displacements  other  than  axial; 

ii.  analyses  of  the  spine  without  head  or  hips; 

iii.  corresponding  analyses  on  rasad-King  (1974)  model. 

The  last  analyses  was  performed  to  insure  that  the  low 

frequency  response  was  not  an  idiosyncracy  of  our  model.  The 
Prasad-King  model  is  somewhat  stiffer,  for  evidently  the  preload 
is  included  in  the  material  properties.  However,  the  differences 
are  quantitative  rather  than  qualitative:  the  Prasad-King  model 
low  frequency  content  is  also  entirely  flexural,  and  the  axial 
mode  frequency  is  even  higher.  It  becomes  clear  from  these  results 
that  contrary  to  widespread  notions,  the  axial  mode  is  never  near 
10  cps,  and  in  fact  hand  calculations  show  that  such  axial 
frequencies  are  almost  impossible.  Evidently,  the  peaks  in  the 
impedance  curves  found  in  axial  harmonic  oscillations  of  the  human 
body  are  due  to  the  parametric  excitation  of  bending  modes  ^n  the 
spii'ic.  This  liypothesis  will  be  explored  further. 


4.  Cervical  Spine  Model 


The  simulations  reported  in  this  section  were  conducted  with 
the  detailed  model  of  the  cervical  spine  described  in  Chapter  IV. 
Additional  beam  elements  which  act  only  in  bending  and  interconnect 
the  primary  nodes  were  included  in  the  cervical  spine.  These  were 
added  to  represent  the  stiffness  of  the  cervical  musculature; 
//ithout  them  it  was  found  impossible  to  maintain  stability  of  the 
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cervical  spine. 

Because  of  the  high  stiffness  of  the  facet  models,  including 

the  cervical  region  in  the  model  reduces  the  stability  limit  of 

the  explicit  integration  procedure  by  an  order  of  magnitude. 

Explicit  integration  would  therefore  require  8000  time  steps  for 

an  80  msec  simulation.  Because  of  the  small  bandwidth  of  this 

model,  implicit  integration  was  quite  suitable.  A time  step  of 
-4 

5x10  seconds  was  used,  and  160  time  steps  were  needed  for  a 80 
msec  simulation.  Although  the  implicit  procedure  increases  the 
stability  limit,  .it  is  not  well-suited  for  contact  problems,  such 
as  the  seatback  interaction.  Thus  a modified  seatback,  which 
restricted  both  posterior  and  anterior  sagittal  plane  motion  was 
used.  In  all  cases  ejection  acceleration  was  prescribed  as  lOG^, 
with  an  onset  rate  of  714  g/sec  over  the  first  14  msec  to  the 
maximum  acceleration. 


Symmetric  Head  Loading 

Comparison  of  Detailed  and  Simple  Cervical  Spine  Models.  The 
simple  model  of  the  cervical  spine  consists  of  a single  beam 
element  which  connects  the  top  of  T1  to  the  head  mass.  In  both 
the  detailed  and  simple  models  the  head  mass  is  the  same  but  the 
detailed  model  includes  the  inertia  of  the  cervical  spine  at  each 
vertebral  level.  A comparison  of  maximum  force  levels  in  the 
upper  thoracic  region  predicted  by  the  two  models  shows  an  in- 
crease in  axial  forces  and  a reduction  in  sagittal  moments  of  the 
intervertebral  disc  for  the  detailed  model.  The  increase  in 
axial  loads  can  be  attributed  to  the  additional  inertia  load  of 
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the  cervical  vertebrae,  however  the  change  in  sagittal  moments 
could  not  be  explained.  This  reduction  in  bending  response  also 
results  in  an  overall  reduction  of  articular  facet  forces  in  the 
detailed  model. 

Compressive  axial  forces  in  the  cervical  region  ranged  from 
1900  Newtons  at  C7  to  1000  Newtons  at  C2  whereas  the  single  beam 
element  in  the  simple  model  has  a maximiam  of  1000  Newtons  since 
only  the  inertia  of  the  head  was  included.  Sagittal  plane  moments 
in  the  cervical  region  are  distributed  between  the  intervertebral 
discs  and  the  added  beams  which  model  the  musculature  with  the 
added  beams  c enerally  having  twice  the  disc  moment.  Both  the 
disc  aj.d  added  beams  have  comparatively  large  moments  at  the 
bottom  and  top  of  the  cervical  region,  while  the  largest  facet 
forces  occur  in  the  center. 

Increased  Head  Loading.  This  study  was  made  to  evaluate  the 
response  of  the  cervical  spine  when  a 0.9  kg  mass  is  added  to  the 
center  of  gravity  of  the  head  mass.  The  cervical  force  levels 
with  the  added  mass  are  significantly  greater.  Axial  disc  forces 
increased  14%  at  C2  to  6.7%  at  Tl  and  sagittal  plane  moments 
averaged  a 6%  increase  with  a corresponding  average  increase  of 
10%  in  articular  facet  forces.  Hence  increasing  the  inertial 
load  of  the  head  by  14%  results  in  an  almost  equal  redistribution 
of  the  added  load  in  the  axial  component  of  the  disc  and  the 
articular  facets. 

In  both  of  these  simulations  the  head  mass  displaced  about 
0.5  cm  in  the  posterior  direction  as  a result  of  locating  the 
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Table  17.  Cervical  Forces  with  Increased  Head  Mass 


center  of  mass  slightly  posterior  to  the  line  of  action  of  the 
cervical  elements.  This  position  of  the  head  mass  was  chosen  to 
insure  that  the  head  would  contact  the  seatback.  An  anterior 
motion  of  the  head  can  be  simulated  by  shifting  the  location  of 
the  head  mass  forward.  However  this  would  require  a provision  for 
possible  chin-chest  contact. 


Ecc.iUtric  Head  Loading.  In  this  study,  the  0.9  kg  mass  was  placed 
10.16  cm  directly  to  the  left  of  the  head  mass  center  of  gravity. 
This  results  in  a 1.41  cm  shift  in  the  center  of  gravity  of  the 
entire  head  mass  to  the  left  of  the  sagittal  plane.  Comparing 
the  results  for  eccentric  and  symmetric  added  mass  shows  almost 
no  change  in  the  intervertebral  disc  forces  with  the  exception  of 
frontal  plane  moments.  However,  as  expected  significant  frontal 
plane  moments  appear  when  the  head  mass  is  eccentric. 

Of  particular  interest  are  the  changes  in  the  articular 
facet  load  distributions.  The  largest  changes  occur  in  the  upper 
three  facet  levels,  C3-C2,  C4-C3,  and  C5-C4,  where  the  facet 
forces  on  the  same  side  as  the  eccentric  mass  are  increased  an  aver- 
age of  37%,  while  the  right  side  facet  forces  are  decreased  by 
about  20%.  The  lower  two  facet  pairs,  C6-C5  and  C7-C6,  exhibit 
a reversal  in  load  distribution  with  the  left  side  facet  loads 
decreasing  by  11%  and  the  right  side  facet  loads  increasing  by 
about  5%.  This  redistribution  of  facet  loads  combined  with  the 
frontal  plane  moments  in  the  intervertebral  disc  enable  the  cervi- 
cal spine  to  carry  almost  all  of  the  moment  due  tc  the  eccentric 
head  mass.  Time  histories  comparing  the  eccentric  and  symmetric 
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case  are  shown  in  Figs.  38  to  42. 

In  the  above  results,  the  added  beams  had  no  frontal  plane 
stiffness.  When  frontal  plane  bending  stiffness  is  included  in 
these  beams,  part  of  the  frontal  plane  moment  due  to  the  eccentric 
head  mass  is  carried  in  these  elements.  In  this  case  all  of  the 
facet  loads  on  the  side  of  the  eccentric  mass  are  increased  by 
about  20%.  The  right  side  facet  loads  were  decreased  from  24%  at 
C3-C2  to  no  change  at  C7-C6.  Generally  by  including  the  frontal 
plane  bending  in  the  additional  beams  both  the  disc  moments  and 
facet  loads  are  decreased. 

A comparison  of  the  simple  cervical  spine  model  with  the 
detailed  model  shows  the  lateral  displacement  of  the  head  is 
0.44  cm  for  the  simple  model  and  0.52  cm  for  the  detailed  model 
at  80  msec.  Also  the  frontal  plane  moment  in  the  simple  model  is 
an  order  of  magnitude  greater  than  frontal  plane  moments  predicted 
by  the  detailed  model.  Thus  ascribing  all  of  the  frontal  plane 
moment  to  the  articular  facets,  as  was  done  in  the  first  section 
of  results  for  eccentric  head  loading,  leads  to  an  exaggerated 
estimate  of  the  facet  loads.  Maximum  head  accelerations  were  18.9 
for  the  simple  model  and  20. 4G  for  the  detailed  model  which  are 
similar  to  the  results  for  the  symmetric  head  loading.  Frontal 
and  lateral  views  of  the  eccentrically  loaded  detailed  model  at 
80  msec  are  shown  in  Figs.  43  and  44. 
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MOMENT  - NEWTON-METERS 


C5-C4  SAGITTAL 


PLANE  MOMENT 


TIME  - SECONDS 

Figure  40.  Comparison  of  C5/C4  sagittal  plane  moment  for 
symmetric  and  eccentric  head  mass. 
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APPENDIX  I 

COMPUTER  PROGRAM  DESCRIPTION 

1.  Introduction 


The  program  package  developed  here  consists  of  two  distinct 
programs;  an  analysis  program  for  predicting  the  dynamic  response 
of  the  human  body  under  prescribed  loads  or  accelerations  and  a 
graphics  package  for  depicting  deformed  and  undeformed  anatomical 
configurations . 

The  techniques  employed  in  the  analysis  program  have  been 
described  in  Chapters  2 and  3.  In  the  following  Sections,  the  in- 
put formats  for  this  program  and  other  information  needed  for  its 
use  are  given.  The  program,  in  addition  to  standard  printer  out- 
put, has  provisions  for  Calcomp  and  printer  plot  graphical  output 
of  time  histories  of  responses  such  as  displacements,  velocities, 
accelerations,  forces,  stresses  and  strains. 

However,  the  graphical  depictions  of  the  anatomy,  such  as 
Figures  6,  7,  8,  and  9 are  obtained  by  a separate  program,  which 
is  described  in  Sections  12  and  13.  The  input  for  this  graphics 
package  can  be  automatically  generated  by  the  analysis  program  in 
punched  card  form.  In  addition,  the  graphics  program  can  be  used 
independently  for  other  studies. 

Both  programs  were  developed  on  an  IBM  370/158  computer  sys- 
tem. They  are  completely  written  in  FORTRAN  IV,  and  with  small 
changes  can  be  run  on  CDC  and  UNIVAC  computers.  For  the  analysis 
program,  both  the  running  time  and  required  core  storage  are 
strongly  problem  dependent.  A minimum  of  about  320  k bytes  (80,000 
words)  is  nppdpd  for  small  problems,  while  models  which  include 
the  rib  cage  and  certain  implicit  models  require  about  520  k bytes 
(140,000  words).  Running  time  is  3 to  5 minutes  for  the  simple 
models,  10  to  20  minutes  for  the  complex  models  on  the  IBM  370/158. 
For  the  graphics  package,  about  320  k bytes  (80,000  words)  of  core 
storage  are  needed;  running  time  is  about  1 to  2 minutes. 

The  programs  were  developed  as  a research  tool.  An  effort 
has  been  made  to  maintain  modularity  of  subroutine  functions  so 
that  additional  features  may  be  added.  Care  has  been  taken  to  in- 
sure maximum  versatility  and  usability,  but  because  of  the 
evolutionary  character  of  the  development,  these  programs  do  not 
have  the  ease  of  input,  extensive  internal  error  checks,  general- 
ized applicability  of  user-oriented,  general  purpose  program.s. 
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2 , Input 

, Data 

Card  1: 

TITLE  CARD 

(20a4) 

Cols. 

FORTRAN 

Name 

Description 

1-80 

TITLE 

Any  80  alphanumeric  characters  to  identify 
the  problem;  these  characters  will  be 
printed  as  a heading  to  the  output. 

Card  2 ; 

PARAMETER  CARD  ( 81 5 , ElO . 6 , 15) 

1-5 

NNODE 

Number  of  nodes  in  the  model  (includes 
primary  and  secondary  nodes,  but  does  not 
include  orientation  nodes) . 

6-10 

NPRI 

Number  of  primary  nodes  in  the  model. 

11-15 

NAXOR 

Number  of  axis  orientation  nodes  in  the 
model;  these  are  used  to  determine  the 
orientation  of  the  local  y axis  of  the 
beam  elements  (see  Section  5) . 

16-20 

NELE 

Number  of  elements  in  the  model. 

21-25 

NUMMAT 

Number  of  different  element  section  and 
material  types;  each  group  of  Card  4 con- 
stitutes a section-material  type. 

26-30 

NUMDIS 

Number  of  nodes  at  which  any  displacement 
components  are  specified  either  a zero  or 
nonzero  value. 

31-35 

MXSTEP 

Number  of  time  steps  to  be  taken. 

36-40 

NDGREE 

Number  of  degrees  of  freedom  per  node 
(should  De  six) . 

41-50 

DELT 

Time  increment. 

51-55 

NODMAX 

Largest  node  number  in  model  (default 
value:  NODMAX  = NNODE  + NAXOR). 

Card  3;  PROGRAM  CONTROL  CARD  (1615) 

1-5  KONTRL(l)  Global  or  local  coordinate  option  (see 

Section  9) . 

KONTRL(l)  = 0:  All  nodes  are  input  in 
global  coordinates. 
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Cols. 


FORTRAN  Namfc 


Description 


KONTRL(l)  = 1:  Secondary  nodes  are  input 
in  the  local  coordinates  of  the  associated 
primary  node  (rigid  linkage) . 

6-10  KONTRL{2)  Print  option. 

K0NTRL{2)  = 0:  Print  time  h stories  of 
stress,  displacement,  etc. 

KONTRL{2)  = 1:  Omit  printing  of  time 
histories. 


11-15  KONTRLO) 


16-20  KONTRL(4) 

21-25  K0NTRL(5) 

26-30  KONTRL(6) 


Control  parameter  for  initial  body  axes. 

K0NTRL{3)  = 0:  Initial  orientations  (time = 0) 
of  body  coordinates  bj^  of  nodal  masses  are 
taken  to  be  the  principal  axes  as  found  by 
an  eigenvalue  routine. 

KONTRL{3)  = 1:  Initial  orientations  of  body 
coordinates  are  coincident  with  the 
global  coordinates. 

Number  of  sliding  interface  p.lanes  (see 
Section  3)  . 

KONTRL (5) /lOOO  is  the  beta  parameter  in 
the  Newmark  6 integration  (should  be  zero 
for  explicit  tnte.:(ration  and  250  for  im- 
plicit integration) . 

'IMPLICIT  ONLY)  If  the  number  of  iterations 
for  the  last  step  is  i KONTRL (6)  then  a 
nev/  stiffness  matrix  is  formulated  for  the 
next  step,  otherwise  the  previous  stiff- 
ness matrix  is  used. 


31-35 

KONTRL (7) 

(IMPLICIT 
per  step. 

36-40 

KONTRL ( 8 ) 

(IMPLICIT 

criteria. 

41-45 

KONTRL (9) 

(IMPLICIT 

motion. 

ONLY)  Maximum  number  of  iterations 
ONLY)  EPSLON*1000,  energy  error 
ONLY)  Rotational  equations  of 


KONTRL (9)  = 0;  Angular  velocity  product 
terms  omitted. 


46-50 

KONTRL (10) 

KONTRL (9)  = 1:  Angular  velocity  product 
terms  included. 

(IMPLICIT  ONLY)  Geome.tric  stiffness  matrix. 

KONTRL (10)  = 0;  Omitted. 
KONTRL (10)  = 1:  Included. 

Cols.  FORTRAN  Name 


Description 


51-55 


56-60 


61-65 


66-70 

71-75 


76-80 


KONTRL(ll) 


KONTRL(12) 


KONTRL  (13) 


KONTRL (14) 
KONTRL (15) 


KONTRL (16) 


Secondary-Primary  node  identifier  (see 
Card  5 NODAL  DATA  CARD) . 


KONTRL (11)  = 0:  Default  node  type  is 
secondary  node . 

KONTRL(ll)  = 1:  Default  node  type  is 
primary  node. 

(IMPLICIT  ONLY)  Modal  analysis  option. 

KONTRL(12)  = 0:  No  modal  analysis. 

KONTRL (12)  = 1:  Modal  analysis. 

KONTRL (12)  = 2:  Modal  analysis  with  punched 
card  output  for  plotting  mode  shapes  of 
spine. 

Damping  option. 

KONTRL (13)  = 0:  Critical  damping 
KONTRL (13)  = 1:  Viscous  damping 
Not  used. 


Program  restart  option. 

KONTRL (15)  = 0:  No  action. 

KONTRL(15)  > 0:  Step  at  which  all  informa- 
tion is  written  on  unit  13  for  later  re- 
start of  problem. 

KONTRL (15)  < 0;  Number  of  seconds  left  in 
job  when  all  information  is  written  on 
unit  13. 


Debugging  print  option. 

KONTRI.  (16)  = 0:  No  action. 

KONTRL (16)  > 0:  Step  at  which  current  values 
in  all  arrays  are  printed  out  and  execution 
continues . 


i 


Card  4;  MATERIAL  PROPERTY  CARDS  (15 , / , 6E10 . 4 , / , 6E10 . 4 ) 

Three  cards  are  required  per  material  section.  The  input  format 
depends  on  the  element  for  which  the  data  is  needed,  so  choose 
appropriately  from  A,  B,  C,  or  D. 

Card  4A.  Material-section  property  cards  for  3-D  AXIAL  SPRING  ELEMENTS. 
Card  4A. 1 


Material  type  number;  the  following  material 
properties  apply  to  all  elements,  JE,  with 
NODE(9,JE)  = MTYP.  (See  Card  6:  Element 
cards. ) 


1-5 


MTYP 


Cols. 


FORTRAN  Name 


I i 


it  s 

•'ll 


‘iri 

i f 

\ ^ 


i I 

‘i 


Card  4A.2 

1-10  E(1,MTYP) 


11-20 


21-30 

31-40 


E(2,MTYP) 


E(3,MTYP) 
E(4 ,MTYP) 


Description 


Tension  or  compression  cutoff  option. 

E(1,MTYP)  < 0:  Element  has  stiffness  only 
in  compression. 

E(1,MTYP)  = 0:  Element  has  stiffness  in 
both  tension  and  compressio  i. 

E(1,MTYP)  > 0:  Element  has  stiffness  only 
in  tension. 

Either  Young’s  Modulus  or  axial  stiffness 
can  be  specified,  see  E(7,MTYP)  for 
implementing  this  option.  (See  Section 
8 for  consistent  units.) 

Not  used. 

Slack  in  spring,  expressed  as  a strain 
offset  (set  equal  to  zero  if  no  slack  is 
desired) . 


41-50 

E(5,MTYP) 

Not  used. 

51-60 

E (6, MTYP) 

Not  used. 

Card  4A. 

3 

1-10 

E(7,MTYP) 

Cross  sectional  area.  Note:  If  area  is 
specified  as  negative,  then  E (2, MTYP)  is 
the  axial  stiffness  of  the  spring,  i.e. 

E (2, MTYP)  = AE/L. 

11-20 

E (8, MTYP) 

Not  used. 

21-30 

E (9, MTYP) 

Not  used. 

31-40 

E(10,MTYP) 

Not  used. 

41-50 

E (11, MTYP) 

Damping  factor  in  fraction  of  critical 
damping  or  fraction  of  viscous  damping  for 
these  elements,  (see  K0NTRL(13)) 

51-60 

E(12,MTYP) 

Not  used. 

Card  4B: 

Material-section 

property  cards  for  3-D  LINEAR,  ELASTIC 

RECTANGULAR  BEAM 

ELEMENT 

Card  4B. 

1 

1-5 

MTYP 

Material  type  number;  the  following 

JE,  with  NODE(9,JE)  = MTYP  (see  Card  6, 
Element  Cards) . 
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Cols. 

FORTRAN  Name 

Description 

Card  4B.2 

1-10 

E(1,MTYP) 

Density  (see  Section  8 for  consistent 
units ) . 

11-20 

E (2, MTYP) 

Elastic  modulus. 

21-30 

E(3,MTYP) 

Not  used. 

31-40 

E(4,MTYP) 

Not  used. 

41-50 

E (5, MTYP) 

Local  z cross  sectional  dimension. 

51-60 

E(6,MTYP) 

Poisson's  ratio. 

Card  4B.3 

1-10 

E (7, MTYP) 

Cross  sectional  area. 

11-20 

E (8, MTYP) 

Local  y cross  sectional  dimension. 

21-30 

E (9, MTYP) 

Shear  deformation  parameter  in  the  local 
y direction;  set  to  zero  for  no  shear 
effects  (see  Chapter  I). 

31-40 

E (10, MTYP) 

Shear  deformation  parameter  in  the  local 
z direction;  set  to  zero  for  no  shear 
effects  (see  Chapter  I) . 

'1-50 

E (11, MTYP) 

Axial  damping  factor  in  fraction  of 
critical  damping  or  fraction  of  viscous 
damping  for  these  elements,  (see  K0NTRL(13)) 

51-60 

E (12, MTYP) 

Bending  damping  factor  in  fraction  of 

Card  4C;  ^ 
I 

Card  4C.1 

I- 5 

Card  4C.2 
1-10 

II- 20 


critical  damping  or  fraction  of  viscous 
damping  for  these  elements,  (see  K0NTRL(13)) 

Material- sect ion  'iroperty  cards  for  3-D  SPINAL  DISK  BEAM 
ELEMENT 


E(1,MTYP) 
E (2,MTYP) 


Material  type  number,  the  following 
material  properties  apply  to  all  elements, 
JE,  with  NODE(9,JE)  = MTYP  (see  Card  6 
Element  Cards) . 


Axial  stiffness  (see  Section  8 for  con- 
sistent units). 

Torsional  stiffness. 


( 
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Cols . 


FORTRAN  Name  Description 


21-30 

E (3,MTYP) 

Bending  stiffness  about  the  local  y axis. 

31-40 

E (4,MTYP) 

Bending  stiffness  about  the  local  z axis. 

41-50 

E (5,MTYP) 

Cubic  bending  stiffness  about  local  y axis 

51-60 

E (6,MTYP) 

Cubic  bending  stiffness  about  local  z axis 

Card  4C. 

2 

1-10 

E (7,MTYP) 

Not  used. 

11-20 

E (8,NTYP) 

Not  used. 

21-30 

E (9,MTYP) 

Shear  deformation  parameter  in  the  local 
y direction:  set  to  zero  for  no  shear 
effects  (see  Chapter  I) . 

31-40 

E (10,MTYP) 

Shear  deformation  parameter  in  the  local 
z direction;  set  to  zero  for  no  shear 
effects  (see  Chapter  I) . 

41-50 

E (11,MTYP) 

Axial  damping  factor  in  fraction  of 
critical  damping  or  fraction  of  viscous 
damping  for  these  elements  (see  KONTRL (13) ) 

51-60 

E (12,MTYP) 

Bending  damping  factor  in  fraction  of 
critical  or  fraction  of  viscous  damping 
for  these  elements  (see  KONTRL (13)). 

Card  4D: 

Mate rial- sect ion 

property  cards  for  PRESSURE  VOLUME 

PENTAHEDRON  ELEMENT 

Card  4D. 

1 

1-5 

MTYP 

Ha’.ej.icil  type  number,  the  following 
mater  1 properties  apply  to  all  elements, 
JE,  with  NODE(9,JE)  = MTYP  (see  Card  6, 
Element  Cards) . 

Card  4D- 

O 

4^ 

1-10 

E (1,MTYP) 

Bullc  modulus  (see  Section  8 for  consis- 
tent units) . 

11-20 

E (2, MTYP) 

Not  used. 

21-30 

E (3, MTYP) 

Not  used. 

31-40 

E (4, MTYP) 

Not  used. 

41-50 

E (5, MTYP) 

Net  used. 

51-60 

E (6, MTYP) 

Not  used. 
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Cols. 

FORTRAN  Name 

Description 

Card  4d. 

2 

1-10 

E(7,MTYP) 

Not  used. 

11-20 

E(8,MTYP) 

Not  used. 

21-30 

L(9,MTYP) 

Not  used. 

31-40 

E(10,MTYP) 

Not  used. 

41-50 

E(11,MTYP) 

Not  used. 

51-60 

E(12,MTYP) 

Damping  factor  in  fraction 
damping  for  these  elements. 

Card  5 : 

NODAL  DATA  CARDS 

(I5,4x,Al,7El0.4) 

NNODE  and  NAXOR  cards  are  required;  the  orientation  nodes  must 
follow  all  regular  nodes. 

1-5  N Node  number. 

6-10  NODTYP  Secondary-Primary  node  identifier. 

NODTYP  - S:  specifies  this  node  is 
secondary. 

NODTYP  = P;  Specifies  this  node  is  pri- 
mary . 

(NOTE:  This  method  of  node  type  identifica- 
tion is  used  along  with  KONTRL(ll)  to 
identify  primary  nodes  whose  mass  is  cal- 
culated by  the  program,  i.e.  for  3-D 


rectangular  beams.) 

11-20 

XC  (N) 

X - coordiante 

21-30 

YC(N) 

Y - coordinate 

31-40 

ZC(N) 

Z - coordinate 

41-50 

TMASS (1) 

Translational  mass. 

51-60 

TMASS(2) 

Global  X moment  of 

inertia. 

^xx 

61-70 

TMASS (3) 

Global  Y moment  of 

inertia. 

I 

yy 

71-30 

TMASS (4) 

Global  Z moment  of 

inertia. 

I 

zz 

(See  Section  6 for  description  of  mass  lumping  at  the  nodes.) 


card  6:  ELEMENT  DATA  CARDS  (1515) 

NELF  cards  are  required. 

Card  6A:  Element  Data  Card  for  3-D  AXIAL  SPRING;  3-D  LINEAR,  ELASTIC 
RECTANGULAR  BEAM;  and  3-D  SPINAL  DISK  BEAM  ELEMENTS 

1-5  M Element  number. 
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Cols. 


FORTRAN  Name 


Description 


Node  I 

Node  J (The  local  x axis  is  directed  from 
Node  I to  Node  J) . 

If  node  I is  not  a primary  node,  the  pri- 
mary node  associated  with  node  I (see 
Section  5) • 

If  node  J is  not  a primary  node,  the  pri- 
mary node  associated  with  node  J . 

Not  used. 

Not  used. 

Not  used. 

Node  K,  axis  orientation  node  number  to^ 
be  used  in  orienting  the  element  local  y 
axis;  the  y axis  lies  in  the  plane  of 
nodes  I,  J,  and  K (see  Section  5). 

Material  type  number;  for  each  material 
type,  a set  of  material  cards  must  be  pro- 
vided (see  Cards  4,  Material  Property  Cards). 

Element  type  number;  indicates  whether 
element  M is  a spring,  elastic  beam  or 
disc  beam,  or  pressure  volume  element. 


6-10 

11-15 

NODE (1,M) 
NODE(2,M) 

16-20 

NODE(3,M) 

21-25 

N0DE(4,M) 

26-30 

31-35 

36-40 

41-45 

NODE (5, M) 
NODE (6, M) 
NODE (7, M) 
NODE (8, M) 

46-50 

NODE (9, M) 

51-55 

NODE (10,M) 

56-60 

NODE (11, M) 

61-65 

Monc n M > 

66-70 

NODE (13, M) 

71-75 

NODE (14, M) 

NODE (10, M) 

= 1: 

3-D 

NODE ( 10, M) 

= 2; 

3-D 

NODE (10, M) 

= 3: 

3-D 

NODE (10, M) 

= 4: 

j-D 

Not  used. 

KT^4- 

Not  used. 
Not  used. 


axial  spring, 
linear  elastic  beam. 

spinal  disk  ;eam. 
pressure  volume  element. 


:ard  6B:  Element  Data  Cards  for  PRESSURE  VOLUME  PENTRAHEDRON  ELEMENT 


1-5 

M 

Element 

number 

6-10 

NODE(l,M) 

Generic 

node 

I 

11-15 

NODE (2 -M’ 

Generic 

node 

J 

16-20 

NODE ;3,M) 

Generic 

node 

K 

21-25 

MODE  1 4 , M) 

Generic 

node 

L 

26-30 

NODE (5, M) 

Generic 

node 

M 

31-35 

NODE (6, M) 

Generic 

node 

N 
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Cols . 

FORTRAN  Name 

Description 

36-40 

NODE(7,M) 

Primary  node  associated  with  generic  ' 
nodes  I-K. 

41-45 

N0DE(8,M) 

Primary  node  associated  with  generic 
nodes  L-N. 

46-50 

N0DE(9,M) 

Material  type  number;  for  each  material 
type,  a set  of  material  cards  must  be 
provided  (see  Cards  4,  Material  Property 
Cards) . 

51-55 

NODE (10, M) 

Element-  type  number. 

NODE(]0,M)  = 4:  Pressure  volume  pentra- 
hedron . 

56-60 

NODE ( 11, M) 

Not  used. 

61-65 

NODE (12, M) 

Not  used. 

66-70 

NODE (13, M) 

Not  used. 

71-75 

NODE (14, M) 

Not  used. 

Card  7;  PRESCRIBED  DISPLACEMENT  CARDS  (14 , 611 ,E10. 4) 

NUMDIS  cards;  include  only  if  NUMDIS  > 0. 

1-4  N Node  number  at  which  one  or  more  degrees 

of  freedom  are  specified. 


For  each  degree  of  freedom  of  node  N a value  of  I is  specified,  where 


1 


0 indicates  no  constraint  on  that  degree  of  freedom. 

1 indicates  that  the  displacement  or  rotation  compon- 
ent is  always  zero. 

2 indicates  that  the  displaceiueut  or  rotation  compon- 
ent is  prescribed  in  SUBROUTINE  FREEFD 


For  each  component,  a column  is  provided  as  follows 


5 I 

6 I 


7 I 


8 I 


9 .1 


Refers  to  translational  global  x degree  of 
freedom  of  node  N. 

Refers  to  translational  global  y degree  of 
freedom  of  node  N. 

Refers  to  translational  global  z degree  of 
freedom  of  node  N. 


Refers 

degree 

to 

of 

rotation 

freedom. 

about 

body 

X axis 

Refers 

degree 

to 

of 

rotation 

freedom. 

about 

body 

y axis 

I 


) 


I 

I 
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Cols. 


FORTRAN  Name 


Description 


10  I Refers  to  rotation  about  body  z axis 

degree  of  freedom. 

11-20  ANGLE  Not  used. 

Card  7A:  MODAL  ANALYSIS  CARD  (I4,6ll) 

Include  only  if  KONTRL(12)  > 0. 

1-4  NAXISP  Coordinate  axis  for  printer  plots  of  pri- 

mary node  mode  shapes. 

NAXISP  = 1:  Mode  shapes  plotted  vs.  global 
X axis. 

NAXISP  = 2:  Global  y axis. 

NAXISP  = 3:  Global  z axis. 

For  each  of  the  six  degrees  of  freedom  in  the  entire  model  a value 
of  I is  specified,  where 

1=0  indicates  this  D.O.F.  is  to  be  included  in  the 
modal  analysis. 

1=1  indicates  this  D.O.F.  is  to  be  omitted  from  the 
modal  analysis. 

In  addition,  individual  nodal  degrees  of  freedom  may  be  omitted  from 
che  modal  analysis  by  prescribing  that  D.O.F.  to  be  zero  via  the 
prescribed  displacement  cards  (..ee  Card  7)  . 


5 

I 

Refers 

to 

translational  global  x. 

6 

I 

Refers 

to 

translational  global  y. 

7 

I 

Refers 

to 

translational  global  z. 

8 

I 

Refers 

to 

rotation  about  body  x. 

9 

I 

Refers 

to 

rotation  about  body  y. 

10 

I 

Refers 

to 

rotation  about  body  z. 

Card  8 : 

OUTPUT  CONTROL 

CARD  (4110) 

1-10 

NPFREQ 

Frequency  of  output;  whatever  output  is 
desired  will  be  printed  every  NPFREQ  steps 

11-20 

NPRU 

Number 

of 

motion  output  records . 

21-30 

NPRS 

Number 

of 

stress  output  records. 

31-40 

NPIC 

Number 

values 

of 

at 

complete  (all  motion  and  stress 
one  time  step)  output  pictures. 
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Cols. 


FORTRAN  Name 


Description 


Card  9:  MOTION  OUTPUT  CARDS  (110 , lOx, 10A4 ) 
(NPRU  cards;  only  included  if  NPRU  > 0.) 


1-7  UOUT 


Node  number. 


8 J 


9 K 


10  L 


21-60  GLABEL 


Component  number  of  kinematic  variable 
(displacment , velocity  or  acceleration) 
to  be  output. 

J = 1:  Translation  in  global  x direction. 

J = 2:  Translation  in  global  y direction. 

J = 3:  Translation  of  global  z direction. 

J = 4:  Rotation  about  body  x axis. 

J = 5;  Rotation  about  body  y axis. 

J = 6:  Rotation  about  body  z axis. 

Indicates  whether  record  is  displacement, 
velocity,  or  acceleration. 

K = 0:  Displacement  time  history. 

K = 1:  Velocity  time  history. 

K = 2:  Acceleration  time  history. 

Plot  control. 

L = 0:  No  plot  of  time  history. 

L = 1:  Calcomp  plot  of  time  history. 

L = 2;  Printer  plot  of  time  history. 

L = 3 : Both  calcomp  and  printer  plot  of 
time  history. 

L = 4;  Printer  plot  and  punched  cards  of 
time  history. 

Alpnanumeric  information  to  be  printed 
identifying  time  history  plot. 


Card  10;  STRESS  OUTPUT  CARDS  (110 , lOx, 10A4 ) 

(NPRS  cards;  only  included  if  NPRS  > 0.) 

1-7  SOUT  Element  number 

8-9  M Component  number  (see  Section  11  for 

organization  of  STRS  array) . 

10  N N = 0;  No  plot. 

N = 1;  Calcomp  plot  of  time  history. 

N = 2:  Printer  plot  of  time  history. 
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Cols . 


FORTRAN  Name 


Description 

N = 3:  Both  Calcomp  and  printer  plot  of 
time  history. 

N = 4:  Printer  plot  and  punched  cards  of 
time  history. 

21-60  GLABEL  Alphanumeric  information  to  be  printed 

for  identifying  time  history  plot. 

Card  11:  COMPLETE  OUTPUT  PICTURE  CARDS  (2110) 

(NPIC  cards;  only  included  if  NPIC  >0.) 

I- 10  NPOUT  Time  step  at  which  complete  output  picture 

is  desired. 

II- 20  KON  KON  = 1;  Output  displacements  and  unit 

vectors  at  all  nodes  having  mass  (see 
Section  5) . 

KON  = 2;  Output  above  plus  coordinates  of 
deformed  model  (note;  output  for  rotational 
degrees  of  freedom  is  nonsense  for  this 
case) . 

KON  = 3;  Output  above  plus  velocities  and 
accelerations  at  all  nodes  having  mass. 

KON  = 4;  Output  above  plus  all  local  element 
forces . 

KON  - 6 : Output  above  plus  punched  card  out- 
put of  all  deformed  nodes  having  mass  and 
associated  unit  vectors.  (Used  as  input  for 
3-D  plotting  program.) 

KON  < 0;  Punched  card  output  only. 


3.  Ejection  Seat  Geometry  Subroutine 


The  program  provides  the  capability  of  modelling  arbitrary 
ejection  seat  geometries.  Seat  geometry  is  represented  as  a col- 
lection of  piecewise  linear  planes.  Each  plane  is  defined  by 
three  points  (X.,  Xj,  X.)  whose  coordinates  are  given  in  the  global 
system.  The  positive  normal  of  the  plane  is  then  defined  by  the 
right  hand  rule  applied  to  the  sequencing  of  the  points  (X^^,  X^,  Xj)  . 
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In  order  to  prescribe  the  motion  of  the  plane  and  incorporate 
the  restraint  system,  each  plane  is  design? ceC  by  a primary  node. 

The  motion  of  the  primary  node  is  characterized  by  an  acceleration 
vector.  The  magnitude  of  the  acceleration  vector  is  given  by  an 
acceleration  time  history  which  is  input  through  subroutine  ICIF 
(see  Section  4).  The  direction  is  given  by  specifying  the  direc- 
tion cosines  of  the  acceleration  vector.  Each  primary  node  which 
designates  a plane  may  be  associated  with  an  arbitrary  number  of 
secondary  nodes.  These  secondary  nodes  can  then  be  used  as  points 
of  attachment  for  representing  the  restraint  system  with  any  of  the 
deformable  elements  available  in  the  program.  NOTE ; The  primary 
and  secondary  nodes  for  this  system  are  input  through  the  Nodal 
Data  Cards  (see  Card  5) , also  each  primary  node  must  have  a Pre- 
scribed Displacement  Card  (see  Card  7)  with  the  translational 
degrees  of  freedom  indicated  as  prescribed  and  the  rotational  de- 
grees of  freedom  indicated  as  zero.  The  deformable  elements  re- 
presenting the  restraint  system  are  input  through  the  Element  Data 
Cards  (see  Card  6) . 

The  technique  used  in  the  subroutine  is  to  modify  the  equations 
of  motion  for  those  primary  nodes  of  the  model  that  are  in  contact 
with  the  ejection  seat,  primary  nodes  not  in  contact  are  not  effected. 
The  criteria  for  contact  between  a plane  and  a primary  node  of  the 
model  is  that  both  the  relative  displacement  and  acceleration  in 
the  normal  direction  of  the  plane  be  decreasing,  i.e.  the  primary 
node  is  moving  towards  the  plane. 

Input  data  for  subroutine  SLIDER;  this  data  is  placed  after 
Card  11  of  READIN  input  (Note:  data  for  subroutine  ICIF  to  specify 
the  motion  of  the  planes  must  follow  this  data) . For  each  plane 
(I  = 1 to  number  of  planes)  the  following  sequence  of  five  input 
cards  are  required: 

Card  Al:  Plane  Identification  Card  (2I5,6E10.0) 


1 i 

Cols . 

FORTRAN  Name 

Descrii 

ption 

* 1 
i i 

1-5 

NPNO(l) 

Primary  node  number 

designating  the  plane 

1 1 

|i 

6-10 

NASN ( I ) 

Total  number  of  primary  nodes  of  the 
model  which  may  contact  plane  I . 

11-20 

DTCOS(l,I) 

Direction  cosine  of 
with  respect  to  the 

acceleration  vector 
qlobal  x-ayis. 

h 

- 1 

i \ 

21-30 

DIC0S(2,I) 

Direction  cosine  of 
with  respect  to  the 

acceleration  vector 
global  y-axis. 

li 

r ! 

31-40 

DIC0S(3,I) 

Direction  cosine  of 
with  respect  to  the 

acceleration  vector 
global  z-axis. 

41-50 

SEATK(1,I) 

Linear  stiffness  of 

elastic  plane  I. 
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Cols.  FORTRAN  Name  Description 

51-60  SEATK(2,I)  Cubic  stiffness  of  elastic  plane  I. 

61-70  VDAMP(I)  Fraction  of  viscous  damping  for  plane  I. 

Card  A2 ; Contacting  Primary  Node  Numbers  Card  (1015) 

1-5  NCP(J,I)  Node  number  of  primary  node  J which  may 

contact  plane  I;  where  J=1  to  NASN(I). 

NOTE:  If  a consecutive  sequence  of  primary  node  numbers  are  associated 
with  a plane,  the  subroutine  will  generate  the  intermediate  primary 
node  numbers.  The  option  to  generate  intermediate  primary  node  num- 
bers is  indicated  by  specifying  NASN(I)  = 0 on  Card  Al  and  specifying 
the  first  and  last  primary  node  numbers,  NCP(1,I)  and  NCP(2,I),  on 
Card  A2.  This  sequence  must  be  ascending  NCP(1,I)  < NCP(2,I). 

Card  A3;  Plane  Location  Cards  (3E10.0) 

(3  cards)  - for  each  of  the  three  points  X2,  X^)  specify 


1-10 

XI(1) 

Global 

X 

coordinate 

of 

point 

I. 

11-20 

XI(2) 

Global 

Y 

coordinate 

of 

point 

I. 

21-30 

XI(3) 

Global 

Z 

coordinate 

of 

point 

I. 

4.  Subroutine  ICIF;  Cubic  Interpolation  Motion  Record  Subroutine 


The  program  in  its  present  form  can  treat  arbitrary  vertical 
motion  input  of  the  hips  and  seat.  The  motion  input  can  be  speci- 
fied as  either  a displacement,  velocity  or  acceleration  vs.  time 
curve.  For  the  latter  two,  the  program  automatically  integrates 
the  record  once  or  twice,  respectively,  to  obtain  a displacement 
history,  which  is  then  used  to  drive  the  model. 

The  motion  record  is  specified  by  an  arbitrarily  spaced  set 
of  points  (t^,  fj^,  i = 1 to  n,  where 

n - is  the  number  of  points  specified  for  the  motion  record. 

t^  = are  the  times  at  which  the  motion  function  (displacement, 
velocity,  or  acceleration)  is  specified. 

f^  = value  of  the  motion  function  (displacement,  velocity,  or 
acceleration)  at  time  t.. 

f^'  = value  of  the  derivative  of  the  motion  function  with  respect 
to  time  at  time  t^. 

A cubic,  continuously  differentiable  interpolation  function  is 
used  to  approximate  the  motion  record  between  prescribed  points. 
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SUBROUTINE  ICIF  (TIME,  VALUE,  L'il) 

TIME  - Time  at  which  displacement  is  to  be  evaluated. 

VALUE  - Upon  return  contains  value  of  displacement. 

Ml  - Number  of  integrations  to  be  performed, 

NI  = 0,  no  integration,  displacement  data  was  input. 

NI  = 1,  one  integration,  velocity  data  was  input. 

NI  = 2,  two  integrations,  acceleration  data  was  input. 

Input  data  for  subroutine  ICIT;  this  data  is  placed  after  Card  11 
of  RZADIN  input. 


Card  Al:  Initial  Parameter  Card  (I5,5X,2E10. 0) 


Cols , 


11-20 

21-30 


FORTRAN  Name 


NPTS 


Description 


Number  of  points  where  the  motion  functions 
value  and  derivative  are  specified,  n. 

Integration  constant  for  first  integration. 

Integration  constant  for  second  integration. 


Card  A2;  Function  Specification  Cards  (3E10.0) 
(NPTS  Cards)  - for  each  point  i specify 


1-10 

11-20 

21-30 


Time  t . 

1 

Motion  function  value  at  time  t^,  f(t^). 

Derivative  of  function  at  time  t.  , f.'(t.). 

1 1 1 


NOTE:  The  motion  function  must  start  at  time  zero,  t^^  - 0. 
An  example  of  this  data  is  given  below. 

Example:  lOg  Acceleration  Record 


♦2»3 
13  15 


TIME  (miltisec) 
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i 

t . 
1 

f . 
1 

f . ' 
1 

1 

0.0 

O 

• 

o 

7.006x10^ 

2 

13x10"3 

9.108x10^ 

7.006x10^ 

3 

15x10"^ 

9.809x10^ 

0.0 

4. 

1.0 

9.809x10^ 

0.0 

5.  Primary,  Secondary  and  Axis  Orientation  Nodes 


This  program  permits  the  user  to  construct  models  with  or 
without  rigid  linkages  between  nodes.  For  this  purpose,  two 
types  of  nodes  are  used: 

Primary  nodes  - All  degrees  of  freedom  are  associated  with 
primary  nodes.  If  a primary  node  is  associated  with  a 
rigid  linkage,  it  must  be  at  the  mass  center  of  the  rigid 
linkage.  At  most  one  primary  node  is  allowed  per  rigid 
linkage . 

Secondary  nodes  - These  nodes  connect  the  ends  of  deformable 
elements  with  rigid  linkages;  more  than  one  secondary 
node  may  be  associated  with  each  rigid  linkage.  A second- 
ary node  has  no  independent  degrees  of  freedom. 

If  there  is  no  rigid  linkage  at  the  endpoint  of  an  element, 
the  endpoint  is  a primary  node. 

A third  type^of  node,  an  axis  orientation  node,  is  used  to 
orient  the  local  y axis  for  beam  elements. 

Axis  Orientation  node  - The  two  nodes  associated  with  the 
endpoints  of  a beam  element  define  the  local  x axis  for 
the  beam  element.  In  order  to  define  the  local  y axis  for 
the  beeim  element,  a third  node  is  used  to  define  a plane 
containing  both  the  local  x and  y axes.  The  normal  to 
this  plane  is  the  local  z axis. 


6.  Mass  Lumping 

All  masses  are  associated  with  primary  nodes.  The  mass  of  a 
primary  node  may  either  (1)  be  input  directly  through  Cards  5,  (2) 
be  generated  within  the  nrogram  through  lumping  the  masses  of 
elements  or  (3)  by  a combination  of  these  methods  (1)  and  (2). 
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If  method  (!)  is  desired  exclusively,  the  density  of  all 
elements  on  cards  4.1  should  be  input  as  zero. 

Method  (2)  can  be  used  by  omitting  the  input  of  lumped  masses 
and  moments  of  inertia  on  Cards  4 and  inputting  mass  densities  for 
the  elements;  the  program  then  lumps  half  the  mass  and  moments  of 
inertia  of  each  element  at  the  associated  nodes. 

Method  (3)  calls  for  care  on  the  part  of  the  user  when  rigid 
linkages  are  used.  The  program  sums  translational  masses  of  each 
rigid  linkage  as  input  through  Cards  4 and  the  lumped  masses  of 
elements  connected  to  the  primary  nodes.  It  also  simply  sums  the 
moments  of  inertia  of  the  primary  node  and  the  elements  connected 
to  it.  It  does  not  shift  the  canter  of  mass  or  transform  moments 
to  inertia  by  the  parallel  axis  theorem. 

The  analyst  should  also  be  usre  that  all  independent  degrees 
of  freedom  are  associated  with  some  mass.  Any  degree  of  freedom 
not  associated  with  a mass  will  automatically  be  omitted  from  the 
equations  of  motion  and  so  remain  fixed  throughout  the  temporal 
integration. 


7^ Unit  Vectors  Bl,  B2,  and  B3 

The  rotational  equations  of  motion  for  each  node  (Eulor 
equations)  are  formulated  in  ^ coordinate  system  that  rotates  with 
the  node.  This  coordinate  uystem  is  initially  chooser,  to  cc  .iride 
with  the  principal  axes  of  ^ne  moment  of  inertia  tensor  and  re- 
mains so  throughout  the  deformation,  since  it  rotates  with  the 
mass. 

If  the  prinicpal  moments  of  inertia  at  a node  are  about  the 
global  X,  y,  and  z coordinates,  then  initially  the  unit  vectors 
of  the  rotational  coordinate  system  coincide  with  the  global  co- 
ordinates, i.e.  Bl  is  a unit  vector  in  the  global  x direction,  B2 
is  a unit  vector  in  the  global  y direction,  and  B3  is  a unit 
vector  in  the  global  z direction.  Only  the  components  of  two  of 
the  unit  vectors  need  be  known  at  any  time,  i.e.  Bl  and  B3.  The 
third  is  then  given  by  the  cross  product  B2  = B3  x Bl. 
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8.  Consistent  Units 


It  is  recommended  that  the  program  be  used  with  the  following 
system  of  units  (i.e.  all  input  data  should  be  in  these  units): 


Time 

Second 

(sec . ) 

Length 

Centimeter 

(cm. ) 

Mass 

Gram 

(g. ) 

Force 

Dyne 

(g-cm/sec^ ) 

For  the  convenience  of  the  user  the  local  element  forces  and 
moments  are  internally  converted  and  also  output  in  the  english 
units: 

Length  : Inch  (in.) 

Force  : Pound  force (Ibf.) 

However,  any  consistent  system  of  units  may  be  used  for  the 
input  data.  The  resulting  output  will  then  be  in  the  same  system 
(and  the  identifying  units  will  not  be  germaine)  with  the  excep- 
tion of  the  english  units  for  forces  and  moments,  which  are  inter- 
nally converted. 


9 . Global  or  Locai.  cooruxuaLes  for  Se\^ondary  Nodes 

Two  options  are  available  for  the  input  of  the  coordinates 
of  secondary  nodes : 

Global  (KONTRL(l)  = 0) 

All  nodal  coordinates  (primary,  secondary,  and  axis  orientation) 
are  input  in  the  global  x,  y,  and  z coordinates  system. 

Local  (KONTRL(l)  = 1 ) 

The  nodal  coordinates  of  the  secondary  nodes  associated  with 
a rigid  linkage  (primary  node)  are  input  in  local  x,  y,  and 
components  of  the  associated  rigid  linkage.  The  local  axis 
origin  is  the  primary  node  for  that  rigid  linkage. 
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Core  storage  among  the  arrays  is  allocated  internally  within 
the  program.  This  procedure  which  is  called  dynamic  allocation, 
requires  the  user  to  specify  only  the  dimension  of  a single  array. 
All  vectors  and  matrices  needed  for  the  solution  of  the  problem 
are  then  packed  sequentially  into  this  array.  The  size  of  the 
array  is  problem  dependent,  so  that  for  small  problems  little  core 
storage  is  needed.  If  the  dimension  specified  by  the  user  is  not 
sufficiently  large  for  the  problem,  an  error  message  is  printed 
along  with  the  total  amount  of  storage  needed. 

The  following  are  the  problem  dependent  scalar  values  used 
to  allocate  core  storage: 


Values  Obtained  Directly  from  the  Input  Data 


FORTRAN  Name 

IMODAL 

MXSTEP 

NAXOR 

NDGREE 

NELE 

NNODE 

NODMAX 

NOPT 

NPDIC 

NPFREQ 

NPIC 

NPRI 

NPRS 

NPRU 

NUMDIS 

NUMMAT 


Description 

Modal  analysis  option 

Maximum  number  of  time  steps. 

Number  of  axis  orientation  nodes. 

Maximum  number  of  degrees  of  freedom 
per  node. 

Number  of  elements 

Number  of  primary  and  secondary  nodes 
in  the  dataset. 

Maximum  node  number  used  in  the  dataset. 
Number  of  sliding  planes. 

Number  of  uegiees  of  freedom  whose  motion 
is  prescribed  as  nonzero. 

Frequency  of  output;  whatever  output  is 
desired  will  be  printed  every  NPFREQ  steps. 

Number  of  complete  output  pictures. 

Number  of  primary  nodes. 

Number  of  stress  output  records. 

Number  of  motion  output  records. 

Number  of  nodes  at  which  any  displacement 
components  are  specified  either  zero  or 
nonzero . 

Number  of  different  element  section  and 
material  types. 
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Values  Calculated  Internally  from  the  Input  Data 


FORTRAN  Name 

LBKE=NBE* (MUDE+1) 

- ( (MUDE+1) *MUDE)/2 

MEQ=NDGREE*NNODE 


MUD= (MAX(I-J)+1) 
*NDGREE-1 

MUDE=NBE-1 


NBE=NDE*NPRI 


NDE 


NODET=NNODE+NAXOR 


NPLOT=NPRU+NPRS 

NPTS=MXSTEP/ 

NPFREQ+2 

N1=NDGREE*NPRI 


Description 

Required  storage  of  the  stiffness  matrix 
for  modal  analysis. 

Maximum  number  of  degrees  of  freedom  in 
the  mesh. 

Number  of  nonzero  upper  codiagonals  in 
banded  stiffness  matrix. 

Number  of  nonzero  upper  codiagonals  in 
banded  stiffness  matrix  for  modal  analysis. 

Number  of  degrees  of  freedom  in  the  mesh 
for  modal  analysis. 

Number  of  degrees  of  freedom  per  node  for 
modal  analysis  (See  Card  7A) . 

Total  number  of  nodes  in  the  mesh  (primary 
and  secondary  and  axis  orientation) . 

Number  of  time  histories  to  be  plotted. 

Number  of  points  to  be  plotted  for  each 
time  history. 

Number  of  independent  degrees  of  freedom 
in  the  mesh. 


The  following  is  a list  of  the  array  names  and  their  storage 
requirements  which  are  dynamically  allocated: 


FORTRAN  Name 


Size 


Description 


A 

NPTS 

Temporary  storage  of  ordinate  values 
for  plotted  output. 

AL 

NELE 

Element  lengths. 

AO 

N1 

Old  nodal  accelerations. 

AUX 

N1 

Deformed  nodal  coordinates. 

A1 

Nl 

Current  nodal  accelerations. 

BIGK 

Nl* (MUD+1) 
- ( (MUD+1) * 
MUD)/2 

Total  global  stiffness. 

BLAMB 

9*NPRI 

Nodal  body  vector  transformations. 

DICOS 

9*NELE 

Element  vector  transformations. 

DICOSP 

9*N0PT 

Sliding  plane  acceleration  vector 
direction  cosines. 
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FORTRAN  Name 


Size 


.1 


Description 


E 

12*NUMMAT 

Material  and  section  properties. 

EFFLAS 

N1 

Old  effective  force. 

EFFORC 

N1 

Current  effective  force. 

EVAL 

NBE*IMODAL 

Eigenvalues  from  modal  analysis. 

EVEC 

NBE*NBE* 

IMODAL 

Eigenvectors  from  modal  analysis. 

FERROR 

N1 

Error  force. 

FEXOLD 

N1 

Old  external  force. 

FINOLD 

2*N1 

Old  nodal  internal  force. 

FINT 

2*N1 

Nodal  internal  force. 

FORCD 

N1 

Current  external  force. 

GLABEL 

10*NPLOT 

Labels  for  plotted  output. 

INDEX 

NELE+1 

Index  to  the  STRS  array. 

INMESH 

NODMAX 

Internal  node  number  locator. 

IX 

7*NELE 

Element  connectivity. 

MESHIN 

NODET 

Mesh  node  number  locator. 

NASN 

NOPT 

Number  of  primary  nodes  which  may 
contact  the  sliding  plane. 

NCP 

10*NOPT 

Node  numbers  of  primary  nodes  which 
may  contact  the  sliding  plane. 

NODDIS 

NUMDIS+1 

Nodal  fixities. 

NPNO 

NOPT 

Primary  node  number  of  sliding  plane. 

NPOUT 

2*NPIC 

Picture  output. 

NTYPE 

2*NPLOT 

Plot  type  identifier. 

PSU 

NPTS*NPLOT 

Output  values  for  time  histories. 

SEATK 

2*N0PT 

Stiffness  coefficients  for  elastic 
sliding  plane. 

SMASS 

Nl 

Nodal  masses. 

SOUT 

NPRS 

Element  output  code. 

SQMASS 

NBE* IMODAL 

Square  root  of  nodal  masses. 

STFLAS 

LBKE 

Total  global  stiffness  for  modal 
analysis. 

S TOREK 

2*MUD* 

NPDIS 

Storage  of  stiffness  for  prescribed 
displacements . 

STROLD 

(Size  is 

Old  element  information. 

determined 
in  sub- 
routine ASSBLE) 
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FORTRAN  Name 


Size 


Description 

Current  element  information. 


3; 


STRS 

(same  as 
STROLD) 

SS 

NPRS 

T 

NPTS 

TBLAM 

9*NPR1 

UOUT 

NPRU 

UP 

3*N0PT 

UPOLD 

3*N0PT 

UPl 

3*N0PT 

UP2 

3*N0PT 

UU 

NPRU 

VDAMP 

NOPT 

VO 

N1 

VI 

N1 

XC 

NODET 

xo 

MEQ 

XI 

MEQ 

YC 

NODET 

ZC 

NODET 

Output  force  values. 

Time  values  for  plotted  output. 

Temporary  storage  of  body  vector 
transformations . 

Nodal  output  code. 

Displacement  of  sliding  plane. 

Old  displacements  of  sliding  plane. 

Velocity  of  sliding  plane. 

Acceleration  of  sliding  plane. 

Output  kinematic  values. 

Viscous  damping  coefficient  for 
sliding  plane. 

Old  nodal  velocities. 

Current  nodal  velocties. 

X coordinate  of  node  points. 

Old  nodal  displacements. 

Current  nodal  displacements. 

Y coordinate  of  node  points. 

Z coordinate  of  node  points. 


11.  Organization  of  STRS  Array 
3-D  Axial  Spring  Elements 


Component  Number 

Quantity 

1-3 

4-6 

Original  local  coordinates  of  second- 
ary node  I with  respect  to  its  associat 
ed  primary  node,  r“  = secondary-pri- 
mary. 

Body  components  of  r°  (r°=X  r°). 

7-9 

Original  local  coordinates  of  secondary 
node  J with  respect  to  its  associated 
primary  node. 

10-12 

Body  components  of  r“. 

13 

Element  strain 

14 

Element  stress 

15 

Element  axial  foice 

16 

Length  change  in  element 

3-D  Linear,  Elastic 

Rectangular  Beam  Element 

Component  Number 

Quantity 

1-6 

Body  components  of  original  element 
unit  vectors  ej  and  e°  (e‘’=X^ye°) 

7-9 

Original  local  coordinates  of  secondary 
node  with  respect  to  its  associated 
primary  node,  r ° = secondary-primary. 

10-12 

Body  components  of  r°  (r  ° = x!^r ° ) . 

The  above  components  1-12  apply  to 
element  node  I . 

13-24 

Same  as  above  for  element  node  J. 

25 

Element  strain 

26 

Element  stress 

27 

Axial  force 

28 

Local  y shear  force  v 

29 

Local  z shear  force  v^ 

30 

Torsional  moment 

31 

Moment  about  local  y axis  at  node  I 

32 

Moment  about  local  y axis  at  node  J M_„ 

jy 
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33 

34 

35 

36 


42-46 


Moment  about  local  z axis  at  node  I M 
Moment  about  local  z axis  at  node  J M 


Iz 

Jz 


Empty 

Change  in  length 


37 

Rotation 

about 

local 

X 

axis 

0 

X 

38 

Rotation 

about 

local 

y 

axis 

at  node  I 01, 

39 

Rotation 

about 

local 

y 

axis 

at  node  J 0J, 

40 

Rotation 

about 

local 

z 

axis 

at  node  I 01 

41 

Rotation 

about 

local 

z 

axis 

at  node  J 0 J 

Empty 


3-D  Spinal  Disk  Beam  Element 

Component  Number  Quantity 

1-6  Body  components  of  original  element  unit 

vectors  ej  and  e°  {e®  = X^ye®). 


7-9 


10-12 

13-24 

25 

26 

27 

28 

29 

30 

31 

32 

33-34 

35 

36 

37 


Original  local  coordinates  of  secondary 
node  with  respect  to  its  associated  pri- 
mary node,  r®  = secondary-primary. 

Body  components  of  r®{r®  = X*^r®). 

The  above  components  1-12  pertain  to 
element  node  I . 

Same  as  above  for  element  node  J. 

Axial  force 

Local  y shear  force 

Local  z shear  force  V 

z 

Torsional  moment  M^ 

Moment  about  local  y axis  at  node  I Mj^ 

Moment  about  local  y axis  at  node  J M_ 

jy 

Moment  about  local  z axis  at  .node  I M_ 

iz 

Moment  about  local  z axis  at  node;  J M_„ 

Jz 

Empty 

Strain  at  previous  time  step 

Change  in  length 

Rotation  about  local  x axis  0 

X 
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38 

Rotation 

39 

Rotation 

40 

Rotation 

41 

Rotation 

42-46 

Empty 

eibout 

local 

y 

axis 

at 

node 

I 

®iy 

about 

local 

y 

axis 

at 

node 

J 

®Jy 

about 

local 

z 

axis 

at 

node 

I 

about 

local 

z 

axis 

at 

node 

J 

®Jz 

Pressure  Volume  Penetration  Element 


Component  Number 

Quantity 

1-3 

Orignial  local  coordinates  of  secondary 
node  with  respect  to  its  associated  pri- 
mary node  r°  = secondary  - primary. 

4-6 

Body  components  of  r°(r°  = X*^r°).  The 
above  components  1-6  pertain  to  the 
generic  node  1 of  the  element. 

7-36 

Same  as  above  for  the  other  five  generic 
nodes  of  the  element. 

37 

Element  internal  pressure. 

38 

Element  volume. 

39 

Element  area. 

40 

Element  axial  force. 

12.  Graphics  Program  for  Anatomical  Analysis 


This  program  was  developed  to  aid  the  researcher  in  analyzing 
the  simulated  behavior  of  the  human  spine.  The  vertebrae,  rib, 
sternum,  head  and  pelvis  can  be  plotted  with  symbols  that  are  re- 
presentative of  these  anatomical  elements  so  that  the  configura- 
tion of  the  skeletal  system  can  be  visualized.  In  addition, 
options  are  available  for  cylinders  and  additional  geometrical 
figures.  These  bodies  can  be  rotated  and  displaced  to  any  orient- 
ation in  three  dimensional  space  and  this  space  can  then  be 
rotated  relative  to  the  viewing  points.  The  back  and  right  side 
projections  are  then  plotted.  Hidden  portions  of  certain  figures 
are  removed  so  that  resulting  projections  are  not  ambiguous  due 
to  a large  number  of  hidden  lines. 

In  order  to  facilitate  use  of  the  program,  an  overview  of  its 
procedure  is  given  in  the  following  paragraphs.  Further  details 
on  the  algorithms  and  procedures  are  given  in  subsequent  Sections. 
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The  elements,  as  represented  by  the  program,  are  described 
by  a combination  of  points  connected  by  lines  and  by  geometric 
solids,  such  as  spheres  or  cylinders.  The  location  of  the  points 
of  an  element  are  given  in  an  element  coordinate  system  (x,y,z). 
The  origin  of  this  coordinate  system  is  the  base  point  of  the 
element.  This  data  is  called  element  point  data.  In  addition, 
elements  which  include  geometric  solids  require  element  dimension 
data,  such  as  the  radius  and  the  length  of  the  cylinder. 

The  location  and  orientation  of  each  element  is  described  by 
the  coordinates  of  the  base  point  in  the  global  system  (x,y,z) 
and  the  orientation  of  the  element's  local  coordinate  system.  The 
latter  is  specified  by  giving  the  global  components  of  the  element 
coordinate  unit  vectors,  that  is  e^^y,  ®3y»  ®3z' 

these  are  equivalent  to  the  direction  cosines. 

The  plotting  program  thus  requires  two  types  of  data: 

1.  position  and  orientation  data  of  each  element; 

2.  element  description  data  consisting  of  element  point  data 


and/or  element  dimension  data  for  each  element. 

The  first  type  of  data  must  be  input  in  each  run  of  the  pro- 
gram by  means  of  card  input  (Cards  5 in  Section  13) . For  many 
anatomical  elements,  the  element  description  data  is  permanently 
stored  in  a data  bank  and  is  an  essential  part  of  this  program's 
capabilities.  These  elements  are  called  standard  elements  and 
data  for  these  elements  have  been  obtained  by  measurements  of 
skeletal  segments  and  is  reasonably  representative  of  a typical 
anatomy;  this  data  must  be  stored  on  units  20  and  27.  If  foe  any 
reason  the  user  desired  to  modify  this  element  description  data, 
he  must  input  his  own  description  point  data  on  these  data  banks. 
In  addition,  the  user  has  a set  of  commands  available  that  enable 
him  to  plot  simple  elements  such  as  lines,  cylinders  and  spheres: 
these  are  called  user-defined  elements. 


To  aid  in  the  visualization,  the  depiction  may  be  plotted 
from  any  angle  relative  to  the  global  (x,y,z)  coordinates.  This 
option  is  effected  by  the  command  SHFT  (cards  5) , which  implicitly 
constructs  an  (x',y',z')  system,  so  that  the  (x,y,z)  system  is 
shifted  relative  to  (x',y',z')  by  the  Euler  angles  (p,  6,  and  ip. 

The  command  VIEW  (Cards  5)  then  given  the  option  of  plotting  the 
(x',z‘)  or  the  (y',z')  views,  or  both.  If  only  a lateral  view 
and  an  antero-posterior  view  are  desired,  a SHFT  command  is  not 
necessary.  The  coordinates  of  the  physical  plot  are  always  de- 
noted by  X and  Z,  where  X can  correspond  to  either  x'  or  y',  and 
Z always  corresponds  to  z'.  The  X axis  is  placed  along  the  bottom 
of  the  plot,  and  the  Z axis  is  placed  vertically. 

The  dimensions  of  the  cube  in  the  (x',y'z')  system  which  en- 
closes all  elements  in  the  plot  must  be  specified;  anything  out- 
side this  cube  is  not  plotted.  The  physical  dimensions  of  the  plot 
and  location  of  the  origin  of  the  (x',y',z')  coordinate  system  on 


the  plot  are  given  by  Card  4.  The  plot  is  then  automatically 
scaled  so  that  the  specified  space  cube  in  the  (x',y',z')  system 
fits  on  the  specified  size  of  the  physical  plot,  with  the  origin 
of  the  (x',y',z')  system  at  the  point  that  is  specified.  Note 
that  the  physical  dimensions  of  the  plot  are  specified  in  inches, 
whereas  all  data  for  the  depiction  is  given  in  cm. 


13.  Input 


All  data  except  plot  dimensions  are  in  cm.  Plot  dimensions 
are  in  inches. 


Card  1;  Heading  (20A4) 


Cols.  1-80 

TITLE;  80  alphanumeric  characters  to  be 
heading  on  output. 

used  as 

Card  2;  Type 

of  input  (4X,A4) 

Cols.  5-10 

"VECTOR" 

Card  3 ; Space 

in  (x',y'z*)  coordinates  to  be  plotted 

(6F10.0) 

Cols.  1-10 

x' 

min 

11-20 

x' 

max 

21-30 

y ' 

min 

31-40 

y' 

max 

41-50 

7. ' 

min 

51-60 

z' 

Card  4;  Size  of  plots  and  location  of  origin  (6F10.0) 


cols.  1-10  TOTALX  width  of  plot  (inches)  for  plots  of  (x',z') 

view. 

11-20  TOTALY  width  of  plot  (inches)  for  plots  of  (y',z') 
view.  (Only  needed  is  IDN  >1  in  VIEW 
command,  see  Cards  5) . 

21-30  TOTALZ  height  of  plot  (inches) ; will  be  the  same 
for  (x',z')  and  (y',z')  views. 

31-40  XNEG  distance  of  the  origin  of  the  (x',y',z') 
system  from  left  hand  side  of  the  plot 
(inches)  in  (x',z')  plot. 


nniMmmPiipH 


41-50  YNEG 


51-60  ZNEG 


distance  of  origin  of  (x',y',z')  system  from 
left  hand  side  of  the  plot  (inches)  in 
(y',z')  plot. 

distance  of  origin  from  bottom  of  plot 
(inches) . 


COMMAND  cards ; Each  group  of  COMMAND  cards  plots  a certain  type 
of  element  or  executes  a specific  operation  in  the  graphics  pro- 
cedure. For  each  command,  Cards  5A,  5B,  and  5C  must  be  in 
sequence.  In  some  commands.  Cards  5A  and/or  5C  are  not  needed 
and  should  be  omitted;  this  varies  with  the  command  and  the  user 
should  refer  to  COMMAND  descriptions. 


Cards  5A;  Command  cards  - choose  any  of  the  commands  listed  under 
COMMAND  DESCRIPTIONS,  (212, A4 , 9F8 . 0) 

Cols.  1-2  IDN  Body  number  (only  needed  for  certain 

standard  elements,  see  5-1)  . 

3-4  NP  Number  of  element  data  points;  if  body  is  on 

data  bank  (see  Table  18) , set  NP  = 0 and  the 
element  description  data  need  not  be  input. 
If  an  element  is  not  a standard  element,  set 
NP  = number  of  data  points  required  as 
specified  by  the  element  command  (see  COM- 
MAND description  for  specific  instructions 
for  each  element) . 


5-8 

ID 

Command 

(four  alphanumeric  characters  from 

the  list 

given  in  5 

-1) 

• 

9-16 

P(l) 

Global 

X 

coordinate 

of 

base 

point. 

17-24 

P(2) 

Global 

y 

coordinate 

of 

base 

point. 

25-32 

P(3) 

Global 

z 

coordinate 

of 

base 

point. 

33-40 

XI 

Global 

X 

component 

of 

local 

X unit 

vector , 

e-' 

XX 

41-48 

X2 

Global 

y 

component 

of 

local 

X unit 

A 

vector. 

e^' 

xy 

49-56 

X3 

Global 

z 

component 

of 

local 

X unit 

vector. 

S'- 

xz 

57-64 

Z1 

Global 

X 

component 

of 

local 

z unit 

vector. 

e ^ 
zx 

65-72 

Z2 

Global 

y 

component 

of 

local 

z unit 

vector , 

e^' 

zy 

73-80 

Z3 

Global 

z 

component 

of 

local 

z unit 

vector. 

e^' 

zz 

5B;  Fo] lows 

Cards  5A  for 

certain  commands  as  specified  in 

COMMAND  descriptions. 
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Cards  5C;  (NP/2  cards)  Element  point  coordinate  data.  Needed  only 

for  certain  commands  as  specified  in  COMMAND  description; 
follows  Card  5B;  if  there  is  no  Card  5B  for  that 
COMMAND,  Card  5C  follows  5A. 


Cols. 


X coordinate  of  point  I 
y coordinate  of  point  I 
z coordinate  of  point  I 
X coordinate  of  point  I+l 
y coordinate  of  point  I+l 
z coordinate  of  point  I+l 


1-10 

11-20 

21-30 

31-40 

41-50 

51-60 


COMMAND  Descriptions 


5.1  STANDARD  ELEMENT  PLOT  COMMANDS 


None  of  these  commands  require  a Card  5B  cr  Card  5C. 

HEAD  Plot  the  head.  The  body  number  IDN  is  ignored. 

RIBL  Plot  the  IDNth  left  rib. 

RIBR  Plot  the  IDNth  right  rib. 

STER  Plot  the  sternum.  The  body  number  IDN  is  ignored. 

VERT  Plot  the  INDth  thoracic  vertebra. 

VERL  Plot  the  IDNth  lumbar  vertebra. 

BODY  The  body  with  optional  number  IDN  (see  last  column  in 

Table  18)  will  determine  which  of  the  above  elements  is 
to  be  plotted.  This  option  is  used  only  for  old  data 
sets  which  do  not  include  the  above  names.  BODY  com- 
mands must  precede  all  other  commands. 

5.2  USER-DEFINED  PLOT  COMMANDS 


Cards  5B  and  5C  are  only  needed  as  noted. 


LINC 


LINO 


The  (NP)  points  on  Cards  5C  will  be  plotted  to  form  a 
closed  curve  by  connecting  point  (i)  to  point  (i+l), 
i=l  to  NP,  and  then  connecting  point  (NP)  to  point  (1) . 

The  (NP)  points  on  Cards  5C  will  be  plotted  to  form  an 
open  curve  by  connecting  point  (i)  to  point  (i+l), 
i=l  to  NP. 


LINS 


SPHR 


CONE 


CYLN 


The  NP  points  on  Cards  5C  must  have  positive  x coordinates. 
For  each  point,  an  additional  point  is  generated  by 
mirror  image  in  the  y-2  plane.  The  resulting  set  of 
points  are  then  plotted  as  in  LINC. 

Plot  a sphere  with  optional  line  segment  superimposed. 


Card  5B 
Cols.  1-10 
11-20 
21-30 


x,y,2  coordinates  of  center  of  sphere, 
respectively,  relative  to  base  point, 
usually  0,0,0. 


31-40  x,y,2  coordinates  of  a point  on  the  sphere, 

41-50  respectively,  which  is  used  to  determine 
51-60  the  radius. 


Cards  5C  are  optional  and  can  be  used  to  input  coor- 
dinates x,y,2  or  NP/2  pairs  of  points.  Each  pair  of 
points  is  connected  by  a line  segment. 

A truncated  conic  section  is  to  be  plotted  with 
optional  line  segments  superimposed.  The  base  point 
is  the  center  of  the  bottom  ellipse. 

Card  5B  (6F10.0  Format) 

Cl  - major  half  axis  of  the  bottom  base 
C2  - minor  half  axis  of  the  bottom  base 
C3  - the  height  of  the  body 
C4  - The  major  half  axis  of  the  top  plane 
C5  - the  minor  half  axis  of  the  top  plane 

Cards  5C  are  optional  and  are  used  to  input  line  seg- 
ments just  as  tor  SPHR,  but  only  visible  line  segments 
are  plotted. 

An  elliptical  cylinder  is  to  be  plotted  with  optional 
line  segments  superimposed.  The  base  point  is  the 
center  of  the  bottom  ellipse. 

Card  5B  (3F10.0  Format) 

Cl  - the  major  half  axis  of  the  base 
C2  - the  minor  half  axis  of  the  base 
C3  - the  length  of  the  cylinder 

Cards  5C  are  optional  line  segments  as  in  CONE. 
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5 . 3 COMMANDS  FOR  PROGRAM  EXECUTION  AND  OPTIONS 


(*only  commands  denoted  by  asterisks  must  be  included  for  the  pro- 
gram to  run;  all  others  may  be  omitted  if  the  user  wishes  to  run 
with  the  default  options.) 

GRID  This  command  will  add  a frame  around  the  plot. 

SHFT  This  command  will  rotate  all  bodies  by  Euler  angles 

(in  radians)  (ip,  8,  i|;)  which  are  the  XI,  X2,  X3  values 
on  Cards  5A  respectively. 

VIEW  This  command  determines  the  views  to  be  plotted. 

If  IDN  = 1,  x' ,z'  plane  is  plotted. 

If  IDN  = 2,  y',z'  plane  is  plotted. 

If  IDN  = 3,  both  views  are  plotted. 

If  no  VIEW  card  is  included,  both  views  will  be  plotted. 

DBUG  This  command  will  trace  the  flow  of  the  program  so  error 

conditions  can  be  trac>jd.  This  is  done  by  printing  the 
object  definition  and  tests  performed  by  the  ID/3D  sub- 
routines. This  commcind  will  also  print  the  transformed 
body  point  coordinates. 

MARK  This  command  will  place  a string  of  alphanumeric 

characters  on  the  plot.  The  string  is  taken  from  Card 
5B.  The  orientation  of  the  message  is  determined  from 
the  data  on  Card  5A. 

P(l)  z coordinate  of  the  left  side  of  string. 

P(2)  x'  or  y'  coordinate  of  the  left  side  of  string. 

P(3)  height  of  letters  in  inches. 

XI  angle  at  which  the  message  is  plotted  from  mea- 
sured vertical  axis  in  degrees. 

X2  number  of  characters  in  message. 

X3  viev?  that  message  is  to  be  plotted  on 

X3  = 1 x' ,z ' view 
X3  = 2 y ' , z ' view 

WEDG  This  command  adds  wedging  to  the  vertebrae,  according 

to  the  following; 

IDN  = 1 AP  wedging  only  (AP:  antero-posterior) 

IDN  = 2 AP  and  lateral  wedging 

If  no  WEDG  command  is  read,  only  /ip  wedging  is  included. 

SCAL  This  command  will  plot  a 5 cm  line  on  the  plot  for  in- 

dicating the  scale  of  the  plot. 

NRES  Each  line  segment  is  divided  into  NRES  equal  subdivisions 

to  check  for  hidden  lines.  This  card  may  be  omitted,  in 
this  case  the  value  of  NRES  is  5. 

GO*  This  card  is  used  as  a delimiter  of  the  data  for  each 

problem,  and  must  be  the  last  card  if  more  than  one 
set  of  data  is  to  be  plotted  in  one  run.  A set  of  cards 
for  another  plot  may  follow  this  card.  The  word  GO 
must  be  left  oriented. 
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The  limits  on  the  number  of  bodies  in  the  standard  version 
of  the  program  are: 

maximum  number  of  bodies  MNB  = 60 

maximum  number  of  element  description  data  MWBP  = 1500  x 3 
maximum  number  of  vertebrae  MNV  = 25 
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Table  18 


Standard  Elements 


Body  type 


IDN 

(for  command 
in  left-hand 
column) 


Description 


(Optional  IDN  for 
use  only  when  BODY 
command  is  used.) 


VERT 

1 

Sacrum  (not  plotted) 
1st  thoracic  vertebrae 

1 

2 

VERT 

2 

2nd  thoracic  vertebrae 

3 

VERT 

3 

3rd  thoracic  vertebrae 

4 

VERT 

4 

4th  thoracic  vertebrae 

5 

VERT 

5 

5th  thoracic  vertebrae 

6 

VERT 

6 

6th  thoracic  vertebrae 

7 

VERT 

7 

7th  thoracic  vertebrae 

8 

VERT 

8 

8th  thoracic  vertebrae 

9 

VERT 

9 

9th  thoracic  vertebrae 

10 

VERT 

10 

10th  thoracic  vertebrae 

11 

VERT 

11 

11th  thoracic  vertebrae 

12 

VERT 

12 

12th  thoracic  vertebrae 

13 

VERL 

1 

1st  lumbar  vertebrae 

14 

VERL 

2 

2nd  lumbar  vertebrae 

15 

VERL 

3 

3rd  lumbar  vertebrae 

16 

VERL 

4 

4th  lumbar  vertebrae 

17 

VERL 

5 

5th  lumbar  vertebrae 

18 

RIBR 

1st  right  rib 

19 

RIBL 

1 

1st  left  rib 

20 

RIBR 

2 

2nd  right  rib 

21 

RIBL 

2 

2nd  left  rib 

22 

RIBR 

3 

3rd  right  rib 

23 

RIBL 

3 

3rd  left  rib 

24 

RIBR 

4 

4th  right  rib 

25 

RIBL 

4 

4th  right  rib 

26 
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Table  18  (continued) 


(for  conunand 
in  left-hand 
column) 


Description 


(Optional  IDN  for 
use  only  when  BODY 
command  is  used) 


RIBR 

RIBL 

RIBR 

RIBL 

RIBR 

RIBL 

RIBR 

RIBL 

RIBR 

RIBL 

RIBR 

RIBL 


5th  right  rib 
5th  left  rib 
6th  right  rib 
6th  left  rib 
7th  right  rib 
7th  left  rib 
8th  right  rib 
8th  left  rib 
9th  right  rib 
9th  left  rib 
10th  right  rib 
10th  left  rib 


STER 


HEAD 


sternum 


head 


*IDN  not  needed  for  HEAD  nr  STERN  commands. 


APPENDIX  II.  ELEMENT  STIFFNESS  MATRICES 


In  this  Appendix  the  local  element  stiffness  matrices  and 
corresponding  global  transformation  matrices  are  presented.  The 
elements  are  standard  in  structural  analysis  so  no  derivations  are 
given.  The  local  element  stiffness  consists  of  a tangential  stiff- 
ness plus  a geometric  stiffness 


[^]  = [k.p]  + [kgl  (II. 1) 

which  relate  nodal  deformations  to  nodal  forces  in  the  corotational 
element  coordinates.  The  global  element  stiffness  matrix  is  then 
obtained  by  transforming  the  local  stiffness  by 

[k]  = [T]'^[k]  [T]  (II. 2) 

where  (T]  is  the  matrix  defined  by  Eq.  (2.36).  The  total  global 
stiffness  is  found  by  adding  together  the  global  stiffnesses  of  all 
elements, 


[k] 


e=l 


(II. 3) 


( 0 ) 

where  [il  ] is  the  Boolean  connectivity  matrix  for  the  element.  As 
is  standard  in  finite  element  programs,  the  matrix  operations  in- 
dicated in  Eq.  (II. 3)  are  not  performed  as  matrix  multiplications, 
but  simply  as  additions.  Thus  the  total  global  stiffness  matrix  is 
obtained  by  adding  the  element  global  stiffnesses  into  the  appropriate 
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locations  of  the  total  matrix,  which  of  course  depends  on  the  node 
numbers  of  the  elements,  i.e.  connectivity. 

The  local  element  stiffness  of  spring  elements  is 


1 

0 0 


[k^l 


ki+3k2<S^ 


0 

-1 


0 

0 


sym. 

0 

1 


0 0 0 0 0 
0 0 0 0 0 0 


(II. 4) 


[kcl 


sym. 

0 1 
0 0 0 


0-100 
0 0-10 


1 

0 1 


(II. 5) 


A 

Here  fj^^  is  the  current  axial  force  in  the  element  and  i is  the 
element  length. 

In  the  beam  element,  for  the  purpose  of  saving  space,  it  is  con- 
vfeiiient  to  express  the  tangent  stiffness  matrix  as  a product  of  two 
matrices.  Hence 

Ik^]  = [Q]'^[k*]  [Q]  (II. 6) 
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where 


The  order  of  the  degrees  of  freedom  for  this  stiffness  is  given  by 
Eq.  (2.48). 


Only  the  part  of  the  geometric  stiffness  corresponding  to  the  rota- 
tion of  the  axial  force  is  used,  so  the  geometric  stiffness  of  the 
beam  is  identical  to  that  of  the  spring.  The  hydrodynamic  element 


stiffness  is 


[k^]  = B {EHE}'^ 


(II. 9) 


where  B is  the  tangent  bulk  modulus  and  {E}  is  given  in  Table  1. 

The  geometric  stiffness  is  not  included  for  the  hydrodynamic  element. 

The  transformation  matrices  are  written  as  partitioned  matrices, 
consisting  of  3 x 3 submatrices  defined  by  Eqs.  (2.1)  through  (2.4) 
and  Eq.  (2.18).  The  spring  transformation  is 


T 

[T]" 


[yli 

[0] 

[0] 


[0] 

[0] 


(II. 10) 


The  beam  transformation  is 


[0] 

[0] 

[0] 

[0] 

[0] 

[0] 

[0] 

[0] 

[0] 

[0] 

(II. 11) 


The  hydrodynamic  element  transformation  is  given  by 


